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Approximation  concepts  are  an  effective  approach  for  alleviating  some  of  the 
problems  associated  with  the  direct  use  of  modern  computerized  analysis  techniques  in  an 
optimization  environment.  Recently,  response  surface  approximations  have  gained  popularity 
as  polynomial  approximations  that  are  global  in  nature.  Response  surface  approximations 
shift  the  computational  burden  from  the  optimization  problem  to  the  problem  of 
constructing  the  approximations,  and  accommodate  the  use  of  detailed  analysis  techniques 
without  the  need  of  derivative  information.  Additionally,  response  surface  approximations 
filter  out  numerical  noise  inherent  to  most  numerical  analysis  procedures,  by  providing  a 
smooth  approximate  response  function,  and  simplify  the  integration  of  the  analysis  and  the 
optimization  codes. 

The  present  dissertation  investigates  the  use  of  response  surface  approximations  in 
expensive  structural  optimization  problems  and  aims  to  suggest  techniques  for  improving 


Xll 


both  the  accuracy  of  response  surface  approximations  as  well  as  the  efficiency  with  which  they 
are  constructed.  A stepped  plate  design  problem  is  considered  and  response  surface 
approximations  are  constructed  for  different  failure  mechanisms  using  numerical  experiments 
conducted  with  a finite  element  analysis.  Both  an  isotropic  and  a composite  laminated  plate, 
where  the  change  in  thickness  is  a result  of  internal  ply  drop  off,  are  considered. 

The  proposed  methodology  uses  a combination  of  dimensional  analysis,  higher  order 
response  surface  approximations,  stepwise  regression,  a detailed  error  analysis  and  statistical 
design  of  experiments  to  improve  both  accuracy  and  efficiency.  Dimensional  analysis 
identifies  variables  intrinsic  to  the  problem,  and  thus  reduces  the  number  of  variables  in  the 
resulting  response  surface  approximation.  Stepwise  regression  is  used  to  eliminate 
insignificant  parameters  from  a response  surface  approximation  and  statistical  design  of 
experiments  is  used  to  identify  a small  set  of  data  points  for  constructing  the  response  surface 
approximations,  while  maintaining  accuracy.  The  result  of  the  proposed  methodology  is 
response  surface  approximations  that  are  highly  accurate  over  the  entire  design  space.  The  use 
of  response  surface  approximations  in  expensive  structural  optimization  problems  is 
demonstrated  by  using  the  developed  response  surface  approximations  in  designing  for 
uncertainty.  The  uncertainty  is  modeled  using  fuzzy  set  theory  and  the  optimization  results 
are  presented  in  the  form  of  design  charts. 


xiii 


CHAPTER  1 
INTRODUCTION 


1.1.  Aspects  of  Modern  Engineering  Optimization 

Modern  computerized  analysis  techniques,  in  the  form  of  numerical  simulations, 
enable  the  detailed  analysis  of  complex  systems.  Although  detailed  numerical  simulations 
provide  an  invaluable  tool  in  analyzing  complex  systems,  several  problem  areas  are  associated 
with  the  direct  use  thereof  in  an  optimization  environment.  Most  notable  is  that  the  designer 
is  provided  with  only  a single  design  point  and  lacks  an  overall  perspective  of  the  system 
behavior  within  the  design  space.  Without  this  overall  perspective,  complications  may  arise 
during  the  optimization  procedure. 

Inherent  in  most  numerical  analyses  are  human  errors  resulting  from  dealing  with 
large,  complex  and  cumbersome  codes  (Kaufman  et  al.1  1996),  and  numerical  noise.  Examples 
of  numerical  noise  are  automatic  mesh  generation  procedures  that  provide  finite  element 
meshes  with  an  inconsistent  level  of  fidelity  and  round-off  error  (Burden  and  Faires2  1989, 
p.  10).  Since  human  and  numerical  errors  of  substantial  magnitude  are  usually  detected  and 
corrected,  the  remaining  inconsistencies  have  negligible  influence  on  the  accuracy  of  the 
results,  but  may  cause  non-smooth  functional  behavior  creating  spurious  local  optima,  thus 
making  it  difficult  to  find  the  best  design  (e.g.,  Giunta  et  al.3  1994  and  Dudley  et  al.4  1995). 

Another  problem  area  arises  when  a designer  needs  to  integrate  several  technical 
disciplines,  as  is  typically  required  in  multidisciplinary  design  optimization.  Integrating 
several  software  codes  poses  a significant  organizational  challenge  and  is  a complex  and  often 
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frustrating  task  (Kaufman  et  al.1  1996).  Even  though  the  problem  of  integrating  different 
software  codes  is  amplified  in  multidisciplinary  optimization  applications,  it  may  also  be  an 
important  factor  in  single  discipline  optimization,  where  it  is  difficult  to  link  the  optimization 
algorithm  with  the  analysis  code. 

Finally,  some  optimization  problems  such  as  global  optimization,  multi-objective 
optimization  and  optimization  under  uncertainty  require  a very  large  number  of  analyses  and 
are,  except  for  simple  problems,  often  impractical.  Advances  in  parallel  computing,  which  is 
an  important  aspect  of  modern  computer  systems,  may  help  alleviate  this  problem.  However, 
except  for  some  special  cases  where  multiple  analyses  are  done  simultaneously  with  very  little 
communication  between  the  processors  (coarse  grain  parallelism),  time  consuming  changes  in 
the  analysis  code  are  often  required  in  order  to  take  advantage  of  the  parallelism. 

Detailed  numerical  simulations  thus  exhibit  many  of  the  same  characteristics  of 
physical  experiments.  A physical  experiment  provides  the  designer  with  only  a single  design 
point  in  the  design  space  and  is  prone  to  experimental  error.  Additionally,  it  is  difficult  to 
integrate  physical  experiments  with  an  optimization  algorithm  and  it  is  typically  most  cost 
efficient  to  conduct  physical  experiments  in  large  batches,  consisting  of  a number  of 
independent  experiments  (coarse  grain  parallelism).  These  similarities  between  detailed 
numerical  simulations  and  physical  experiments  are  discussed  in  more  detail  in  a recent  survey 
paper  by  Haftka,  Scott  and  Cruz5  (1998)  that  concentrates  on  optimization  in  relationship  to 
experiments. 

1.2.  Response  Surface  Methodology 

Response  surface  methodology  may  be  summarized  as  a collection  of  statistical  tools 
and  techniques  for  constructing  and  exploring  an  approximate  functional  relationship  between 
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a response  variable  and  a set  of  design  variables.  This  approximate  functional  relationship  is 
typically  constructed  in  the  form  of  a low  order  polynomial,  referred  to  as  a response  surface 
approximation.  Response  surface  methodology  was  originally  developed  for  constructing  and 
exploring  approximate  response  functions  based  on  physical  experiments  and  results  in 
smooth  approximate  response  functions,  thus  effectively  filtering  out  any  noise.  Response 
surface  approximations  allow  easy  integration  between  the  experiments  and  the  optimization 
algorithm,  allow  the  evaluation  of  experiments  in  batches  and  provide  the  designer  with  an 
overall  perspective  of  the  system  behavior  within  the  design  space.  Due  to  the  similarities 
between  detailed  numerical  simulations  and  physical  experiments,  there  has  recently  been 
growing  interest  in  using  response  surface  approximations  within  structural  optimization. 

Response  surface  methodology  includes  methods  for  selecting  data  points  where  the 
experiments  should  be  evaluated  (a  process  known  as  statistical  design  of  experiments), 
methods  for  solving  the  unknown  coefficients  of  a response  surface  approximation  and 
methods  for  evaluating  the  accuracy  of  the  resulting  response  surface  approximation.  The 
unknown  coefficients  of  a response  surface  approximation  are  estimated  from  experimental 
data  points  by  means  of  a process  known  as  linear  regression.  These  coefficients  are  estimated 
in  such  a way  as  to  minimize  the  sum  of  the  squares  of  the  error  terms  (e.g.,  Myers  and 
Montgomery6  1995,  pp.  16-27  and  Box  and  Draper7  1987,  pp.  35-70).  The  accuracy  of  a 
response  surface  approximation  is  expressed  in  terms  of  various  error  terms  and  statistical 
parameters  that  are  representative  of  the  predictive  capabilities  of  the  approximation 
(e.g.,  Myers  and  Montgomery6  1995,  pp.  27-54). 

Statistical  design  of  experiments  is  used  to  select  a small  set  of  data  points  for 
constructing  the  response  surface  approximation  without  compromising  accuracy. 
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Traditionally,  statistical  design  of  experiments  is  performed  using  a minimum-variance 
criterion  that  assumes  all  error  to  be  random  error  or  variance.  For  regularly  shaped  design 
spaces  (i.e.,  rectangular  or  spherical  regions)  there  exist  a number  of  standard  designs  (i.e.,  sets 
of  data  points),  for  example  the  Central  Composite  (Box  and  Wilson8  1951)  and  the  Box- 
Behnken  (Box  and  Behnken9  1960)  designs.  For  an  irregularly  shaped  design  space,  a 
computer-generated  design  must  be  used,  of  which  the  most  popular  employ  the 
Z>Optimality  design  criterion  (e.g.,  Myers  and  Montgomery5  1995,  pp.  364-366). 

Alternatively,  the  data  points  may  be  selected  to  minimize  the  bias  error  of  the 
resulting  response  surface  approximations.  These  designs  are  known  as  minimum-bias  designs. 
Bias  error  occurs  when  a response  surface  approximation  that  is  inadequate  in  modeling  the 
actual  response  function  is  used.  Since  the  exact  functional  form  of  the  response  function  to 
be  approximated  is  rarely  known  and  low  order  polynomials  are  generally  used  as  response 
surface  approximations,  large  bias  error  is  easily  introduced.  For  regularly  shaped  design 
spaces,  analytical  solutions  exist  for  finding  minimum-bias  designs  (Myers  and  Montgomery6 
1995,  pp.  406-421  and  Box  and  Draper7  1987,  pp.  437-442).  Balabanov  et  al.10  (1996) 
successfully  applied  these  analytical  solutions  to  find  a minimum-bias  design  for  a 25- 
dimensional  sphere.  For  irregularly  shaped  design  spaces,  no  closed  form  solution  exists  and 
the  minimum-bias  data  points  must  be  obtained  numerically  (e.g.,  Venter  and  Haftka"  1997). 

An  approach  where  the  response  surface  approximation  is  built  from  data  points 
selected  by  the  optimizer  has  also  recently  been  used  in  structural  optimization.  For  example, 
Harrison,  Le  Riche  and  Haftka12  (1995)  and  Sellar  and  Batill 13  (1996)  used  the  approach  to 
construct  response  surface  approximations  that  are  periodically  updated  throughout  the  design 
process.  The  response  surface  approximation  is  updated  by  adding  new  design  points, 
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generated  by  the  optimizer,  to  the  set  of  design  points  used  for  constructing  the  response 
surface  approximation.  Wujek  et  al.14  (1996)  and  Rodriguez,  Renaud  and  Watson15  (1997) 
applied  the  approach  to  multidisciplinary  design  optimization,  constructing  local  response 
surface  approximations  in  a sequential  manner  based  on  design  points  evaluated  at  the 
discipline  level  of  the  current  design  cycle.  Additionally,  Cramer  and  Carter16  (1997) 
considered  design  points  generated  by  a global  optimization  algorithm.  A statistical  analysis 
of  these  design  points  was  then  used  to  complement  the  global  optimization  process  by  giving 
more  insight  into  the  system  behavior  within  the  design  space.  This  additional  information 
was  then  used  to  influence  the  design  process,  for  example  to  eliminate  design  variables  that 
have  a small  influence  on  the  system  response. 

1.3.  Response  Surface  Approximation  Applications  in  Optimization 

As  discussed  in  the  previous  section,  the  problems  associated  with  the  direct  use  of 
detailed  numerical  simulations  in  an  optimization  environment  may  be  alleviated  by  the  use 
of  approximation  concepts.  When  using  approximation  concepts,  the  computational  burden 
is  shifted  from  the  optimization  problem  to  the  problem  of  constructing  the  approximations. 
A detailed  discussion  of  basic  approximation  concepts  in  structural  optimization  is  presented 
by  Barthelemy  and  Haftka17  (1993). 

It  has  long  been  recognized  in  the  field  of  structural  optimization  that,  compared  to 
some  of  the  direct  methods,  for  example  the  method  of  feasible  directions  (e.g.,  Haftka  and 
Giirdal18  1993,  pp.  182-186),  approximation  concepts  are  an  effective  approach  for  reducing  the 
computational  cost  associated  with  optimization  problems.  Schmit  and  his  co-workers 
(e.g.,  Schmit  and  Farshi19  1974  and  Schmit  and  Miura20  1976)  demonstrated  first  that,  when 
applied  correctly,  approximation  concepts  might  significantly  improve  the  efficiency  of 
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solving  structural  optimization  problems.  The  traditional  approach  is  to  use  local  derivative- 
based  approximations  (typically  some  variation  of  the  Taylor  series  expansion)  that  are 
generally  applied  sequentially,  with  new  approximations  constructed  at  each  new  design  cycle. 
These  local  derivative-based  approximations  are  constructed  from  information  obtained  at  a 
single  design  point,  and  help  to  reduce  the  magnitude  of  the  problem  areas  associated  with  the 
use  of  detailed  numerical  analyses  in  an  optimization  environment.  Unfortunately,  as  a result 
of  their  local  nature  and  sequential  application,  the  problem  areas  are  not  eliminated.  The 
designer  still  lacks  a global  perspective  of  the  response,  human  error  and  numerical  noise  may 
still  cause  the  optimization  algorithm  to  converge  on  spurious  local  optima,  it  is  still  difficult 
to  take  full  advantage  of  parallel  computing  techniques,  and  interaction  between  different 
analysis  codes  is  still  required. 

Multi-point  approximations  are  a variation  of  the  traditional  local  derivative-based 
approximations  that  utilize  both  function  and  derivative  information  from  several  design 
points  in  constructing  an  approximation.  Compared  to  traditional  derivative-based 
approximations,  typically  constructed  from  function  and  derivative  information  obtained 
from  a single  design  point,  multi-point  approximations  are  generally  more  accurate  and  global 
in  nature.  Haftka  et  al.21  (1989)  examined  approximations  based  on  two  and  three  design 
points  and  concluded  that  the  approximations  worked  well  for  interpolation,  while  only 
marginal  improvements  in  accuracy  were  obtained  when  extrapolating.  Fadel  et  al.22  (1990) 
introduced  a two  point  exponential  approximation  that  uses  derivative  information  from  the 
previous  design  point,  while  Snyman  and  Stander23  (1994)  used  only  the  previous  function 
value  in  creating  a two  point  quadratic  approximation.  Additionally,  Wang  and  Grandhi24,25 
(1994  and  1995)  introduced  two  point  approximations  that  use  both  function  and  derivative 
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information  from  a previous  design  point.  Wang,  Grandhi  and  Hopkins26  (1994),  Wang, 
Grandhi  and  Canfield27  (1994)  and  Wang  and  Grandhi28  (1996)  introduced  several  other  multi- 
point approximations,  using  information  from  more  than  two  design  points.  These 
approximations  were  applied  successfully  to  various  truss,  frame,  shape,  and  reliability-based 
optimization  problems.  In  their  work,  the  set  of  design  points  used  for  constructing  the 
approximation  is  expanded  at  the  end  of  each  design  cycle  by  performing  an  accurate  analysis 
of  the  current  optimum  design.  As  the  design  process  proceeds  more  design  points  are  used 
for  constructing  the  approximation,  thus  gradually  improving  the  accuracy  and  expanding  the 
region  of  applicability. 

More  recently,  response  surface  approximations  have  been  used  extensively  for 
obtaining  multi-point  approximations.  Compared  to  traditional  multi-point  approximations 
as  discussed  in  the  previous  paragraph,  the  multi-point  approximations  based  on  response 
surface  methodology  do  not  have  to  interpolate  between  the  design  points  but  use  a least 
squares  procedure.  Multi-point  approximations  based  on  response  surface  methodology  differ 
from  standard  response  surface  approximations,  since  they  are  constructed  from  design  points 
selected  by  the  optimizer.  In  contrast,  the  standard  response  surface  approximations  are 
constructed  from  data  points  obtained  from  statistical  design  of  experiments.  Additionally, 
some  of  these  multi-point  approximations  make  use  of  derivative  information  in  constructing 
the  approximations  (e.g..  Sellar  and  Batill13  1996  and  Rodriguez,  Renaud  and  Watson15,29  1997) 
while  standard  response  surface  approximations  use  only  function  values. 

In  the  literature,  multi-point  approximations  based  on  response  surface  methodology 
are  implemented  in  many  different  formulations.  For  example,  Toropov,  Filatov  and 
Polynkin30  (1993)  used  response  surface  methodology  to  develop  a multi-point  approximation 
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design  methodology.  The  approach  employs  move  limits  during  the  optimization  process  to 
define  a sub-region  about  the  current  optimum  design  point.  A weighted  least  squares 
procedure  is  then  used  to  formulate  an  explicit  local  approximation  near  the  current  design 
point,  according  to  the  move  limit  strategy.  Both  new  design  points  evaluated  within  the 
current  sub-region  as  well  as  several  design  points  evaluated  during  previous  design  cycles  are 
used  to  construct  the  approximations.  The  methodology  is  applied  successfully  to  several 
example  problems,  e.g.,  truss  and  beam  structures  (Toropov,  Filatov  and  Polynkin30  1993), 
performance  optimization  of  small  Stirling  engines  (Toropov,  Markin  and  Carlsen31  1994  and 
Toropov  and  Carlsen32  1994)  and  the  optimization  of  geometrically  nonlinear  shell  structures 
(Van  Keulen,  Toropov  and  Polynkin33  1994). 

Sellar  and  Batill1'  (1996)  also  constructed  multi-point  approximations  based  on 
response  surface  methodology,  but  used  a slightly  different  implementation  specifically 
tailored  to  a multidisciplinary  optimization  environment.  At  each  new  design  cycle,  function 
and  derivative  information  obtained  from  design  points  identified  by  the  optimizer  were 
added  to  the  database  used  for  constructing  the  response  surface  approximation.  The  response 
surface  approximation  is  thus  updated  at  each  new  design  cycle  with  new  information 
obtained  from  design  points  identified  by  the  optimizer.  Harrison,  Le  Riche  and  Haftka12 
(1995)  used  a similar  approach.  They  made  use  of  a genetic  algorithm  to  minimize  the  weight 
of  a hat  stiffened  composite  panel  and  constructed  a response  surface  approximation  based  on 
design  points  identified  by  the  genetic  algorithm.  Their  response  surface  approximation  is 
updated  periodically,  after  a fixed  number  of  iterations  are  completed  by  adding  new  design 
points  identified  by  the  genetic  algorithm  to  the  set  of  design  points  used  for  constructing  the 
response  surface  approximation.  In  their  approach,  the  most  recent  design  points  were  given  a 
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higher  weight  as  compared  to  the  older  design  points,  and  the  most  recent  design  points  thus 
have  a larger  influence  on  the  resulting  response  surface  approximation. 

Venkataraman  and  Haftka34  (1997)  and  Venter,  Haftka  and  Chirehdast35  (1997)  used 
local  response  surface  approximations  constructed  in  a sequential  manner,  in  single  discipline 
design  optimization.  After  each  design  cycle,  a new  response  surface  approximation  is 
constructed  using  statistical  design  of  experiments  to  select  data  points  for  constructing  the 
response  surface  approximation  near  the  current  optimum  design.  Wujek  et  al.14  (1996)  and 
Rodriguez,  Renaud  and  Watson13  (1997)  considered  a similar  approach  in  multidisciplinary 
design  optimization  problems,  but  used  design  points  evaluated  at  the  discipline  level  during 
the  current  design  cycle  for  constructing  the  sequential  response  surface  approximations. 
Rodriguez,  Renaud  and  Watson15  (1997)  made  use  of  trust  region  constraints  to  define  a sub- 
region  near  the  current  optimum  design.  At  each  design  cycle,  a new  response  surface 
approximation  is  constructed  for  the  defined  trust  region.  These  approximations  are 
constructed  so  that  the  function  and  the  first  derivative  match  the  actual  values  at  the  current 
design  point  exactly,  and  thus  only  the  second  order  parameters  are  obtained  using  a least 
squares  procedure.  Rodriguez,  Renaud  and  Watson29  (1997)  constructed  similar 
approximations,  but  used  statistical  design  of  experiments  to  select  the  data  points  for 
constructing  the  response  surface  approximations. 

All  of  the  above  response  surface  approximation  formulations  have  the  advantage  of 
not  requiring  derivative  information  when  constructing  the  approximations.  However,  the 
use  of  derivative  information  is  permitted  and  was  used  in  some  cases  (e.g.,  Ref.  13,  15  and  29). 
Additionally,  these  approximations  have  a tendency  to  eliminate  spurious  local  minima 
resulting  from  human  and  numerical  error,  while  some  of  the  formulations  (Ref.  12  and  13) 
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result  in  approximations  that  are  more  general  in  nature  when  compared  to  the  traditional 
local  derivative-based  approximations.  Although  these  response  surface  approximations  are 
successful  in  alleviating  the  problem  areas  associated  with  using  detailed  numerical  analyses  in 
an  optimization  environment,  the  problem  areas  are  not  completely  eliminated  due  to  the  fact 
that  these  approximations  still  do  not  provide  global  approximations  that  are  valid  over  the 
entire  design  space. 

When  using  polynomial  response  surface  approximations  that  are  global  in  nature, 
however,  the  problems  associated  with  using  modern  analysis  procedures  in  an  optimization 
environment  may  be  eliminated.  Apart  from  providing  the  designer  with  a global  perspective 
of  the  response,  spurious  local  minima  resulting  from  human  and  numerical  errors  are 
eliminated  by  using  the  least  squares  method  to  formulate  smooth  low  order  polynomials. 
These  approximations  are  constructed  independently  and  no  interaction  between  different 
software  codes  is  required.  In  many  cases  the  computational  cost  of  structural  optimization 
problems  may  be  reduced  by  using  response  surface  approximations  that  are  global  in  nature, 
since  the  number  of  function  evaluations  required  to  construct  these  response  surface 
approximations  may  be  less  than  that  needed  to  perform  the  optimization.  Finally,  the  use  of 
parallel  computing  techniques  is  facilitated  in  that  the  response  surface  approximations  are 
constructed  from  independent  data  points,  which  may  be  evaluated  in  parallel  with  very  little 
communication  between  the  different  processors. 

In  recent  years,  many  researchers  have  exploited  the  advantages  of  using  standard 
response  surface  approximations  that  are  more  global  in  nature  to  integrate  numerical  analyses 
in  an  optimization  environment.  With  few  exceptions  (e.g.,  Van  Houten  et  al.36  1995),  almost 
all  the  standard  applications  of  response  surface  approximations  in  structural  optimization 
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employ  quadratic  approximations  which  still  do  not  provide  fully  global  approximations.  For 
example,  Thacker  and  Wu37  (1992)  used  response  surface  approximations  to  reduce  the 
computational  cost  associated  with  reliability-based  optimization  problems  using  probabilistic 
methods,  while  Harrison,  Le  Riche  and  Haftka12  (1995)  reduced  the  cost  of  designing  a 
stiffened  composite  panel  using  a genetic  algorithm.  Van  Houten  et  al.36  (1995),  Giunta  et  al.38 
(1996),  Kaufman  et  al.39  (1996)  and  Venter,  Haftka  and  Chirehdast35  (1997)  showed  that 
response  surface  approximations  are  valuable  in  solving  optimization  problems  with  non- 
smooth functional  behavior,  since  numerical  noise  is  filtered  out  from  the  resulting 
approximate  response  function.  Additionally,  Mistree  et  al.40  (1994)  made  extensive  use  of  the 
global  perspective  of  the  response  over  the  design  space  provided  by  global  response  surface 
approximations,  while  Brown  and  Nachlas41  (1985),  White  et  al.42  (1985)  and  Giunta  et  al.43 
(1995)  used  response  surface  approximations  to  meet  the  organizational  challenge  posed  by 
multidisciplinary  optimization  problems.  More  specifically,  Mason  et  al.44  (1994)  used 
response  surface  approximations  to  integrate  the  analysis  and  optimization  codes  in  a single 
discipline  optimization,  while  Ragon  et  al.45  (1997)  used  response  surface  approximations  as  a 
simple,  yet  flexible,  interface  between  the  global  and  local  design  codes  in  the  design  of  a large 
wing  structure. 

The  Multidisciplinary  Analysis  and  Design  (MAD)  Center  of  Virginia  Tech 
considered  the  aerodynamic-structural  optimization  of  the  High  Speed  Civil  Transport 
(HSCT)  and  made  extensive  use  of  the  fact  that  response  surface  approximations  facilitate  the 
use  of  parallel  computing  techniques  (e.g.,  Giunta  et  al.43  1995,  Balabanov  et  al.10  1996,  Burgee 
et  al.46  1996  and  Kaufman  et  al.39  1996).  They  made  use  of  an  Intel  Paragon  machine  with 
more  than  100  processors  to  perform  a large  number  of  structural  optimizations  and  showed  a 
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substantial  reduction  in  real  time  by  performing  these  optimizations  in  parallel.  The  data 
obtained  from  these  optimizations  were  used  to  construct  response  surface  approximations  for 
the  wing  bending  material  weight  of  the  High  Speed  Civil  Transport. 

Unfortunately,  response  surface  approximations  are  plagued  by  what  is  known  as  the 
curse  of  dimensionality.  The  curse  of  dimensionality  implies  that  the  number  of  data  points 
required  to  construct  a response  surface  approximation  increases  greatly  when  more  design 
variables  are  added  to  the  response  surface  approximation.  A substantial  reduction  in  the  cost 
of  constructing  a response  surface  approximation,  often  equal  to  a few  orders  of  magnitude, 
can  be  realized  by  reducing  the  number  of  design  variables  in  the  response  surface 
approximation.  Identifying  the  variables  that  are  intrinsic  to  the  problem  under 
consideration,  and  using  these  intrinsic  variables  as  design  variables  in  a response  surface 
approximation  is  an  effective  approach  for  reducing  the  number  of  design  variables  in  a 
response  surface  approximation.  The  use  of  intrinsic  variables  has  the  additional  advantage  of 
increasing  the  accuracy  of  the  resulting  response  surface  approximation  (e.g.,  Kaufman  et  al.1 
1996  and  Roux,  Stander  and  Haftka47  1996). 

1.4,  Dimensional  Analysis 

Dimensional  analysis  provides  a methodology  for  identifying  the  intrinsic  variables  of 
a system  by  formulating  the  governing  equations  of  the  system  in  non-dimensional  form. 
This  technique  is  discussed  in  detail  in  Barenblatt48  (1996,  pp.  1-63).  Dimensional  analysis  is 
based  on  the  simple  idea  that  any  physical  law  is  independent  of  the  arbitrary  chosen  units  of 
measurement  used  for  describing  a specific  system.  Instead,  physical  laws  possess  a 
fundamental  property  that  may  be  expressed  in  terms  of  variables  intrinsic  to  the  system.  The 
use  of  dimensional  analysis  as  a general  methodology  for  identifying  the  intrinsic  variables  of  a 
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system  is  justified  by  the  Buckingham  11-Theorem  (Buckingham49  1915),  which  may  be  stated 
as: 

A physical  relationship  between  some  dimensional  (generally  speaking) 
quantity  and  several  dimensional  governing  variables  can  be  rewritten  as  a relationship 
between  some  dimensionless  quantity  and  several  dimensionless  products  of  the 
governing  variables;  the  number  of  dimensionless  products  is  equal  to  the  total 
number  of  governing  variables  minus  the  number  of  governing  variables  with 
independent  dimensions. 

The  governing  variables  (also  referred  to  as  predictor  variables  in  this  dissertation)  refer  to  the 
dimensional  variables  used  in  describing  the  system.  An  independent  dimension  is  defined  as  a 
dimension  that  cannot  be  expressed  in  terms  of  a product  of  powers  of  the  dimensions  of  the 
remaining  variables. 

By  using  the  non-dimensional  intrinsic  variables  obtained  from  the  non-dimensional 
formulation  of  the  governing  equations,  the  number  of  design  variables  in  a response  surface 
approximation  may  be  greatly  reduced.  Additionally,  a non-dimensional  formulation  of  the 
governing  equations  of  a system  results  in  more  general  solutions  that  are  applicable  to  a wider 
range  of  problems.  A good  example  is  the  stress  concentration  design  charts  published  by 
Peterson50  (1953).  These  charts  are  applicable  to  any  isotropic  and  homogeneous  material  and 
are  independent  of  the  absolute  geometry  and  loading  conditions  of  the  structure  under 
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1.5.  Fuzzy  Set  Based  Design  Optimization 
A typical  example  of  an  expensive  structural  optimization  problem  that  will  be 
considered  in  the  present  dissertation  is  designing  for  uncertainty.  Designing  of  uncertainty 
often  occurs  in  the  aircraft  industry,  where  structures  are  designed  that  will  be  built  well  into 
the  future  from  materials  available  then.  Oftentimes,  these  materials  are  not  yet  available, 
forcing  the  designer  to  deal  with  uncertain  material  properties.  For  this  type  of  design 
problem,  little  information  regarding  the  uncertainty  is  known  and  the  uncertainty  is 
typically  based  on  expert  opinion  and  assumptions  made  by  the  designer.  Unlike  probabilistic 
methods,  which  require  a large  amount  of  data  for  modeling  the  uncertainty,  fuzzy  set  theory 
is  capable  of  modeling  uncertainty  using  the  limited  data  available  to  the  designer.  Apart  from 
requiring  a large  amount  of  data,  the  results  obtained  from  probabilistic  methods  are,  in  some 
cases,  very  sensitive  to  the  accuracy  of  the  data  as  well  as  to  assumptions  made  during  the 
modeling  process  (Ben-Haim  and  ElishakofP1  1990,  pp.  11-32).  In  contrast,  fuzzy  set  theory 
caters  to  worst  case  scenarios,  taking  account  of  the  fact  that  the  uncertainty  is  modeled  based 
on  subjective  opinions  and  assumptions  (Maglaras  et  al.52  1997). 

Fuzzy  set  theory  was  introduced  by  Zadeh53  (1965)  as  a mathematical  tool  for  the 
quantitative  modeling  of  uncertainty.  Fuzzy  set  theory  makes  use  of  fuzzy  numbers  to 
represent  uncertain  problem  parameters.  In  order  to  represent  an  uncertain  parameter  as  a 
fuzzy  number,  the  designer  only  needs  to  specify  the  range  of  uncertainty  and  a membership 
function  that  denotes  the  possibility  of  occurrence  of  an  element  in  the  specified  range. 

In  recent  years,  fuzzy  set  theory  has  been  applied  to  a wide  range  of  structural 
optimization  problems.  For  example,  Liu  and  Huang54  (1992)  performed  a fatigue  reliability 
analysis  of  a portal  frame,  Jung  and  Pulmano55  (1996)  considered  the  optimal  plastic  design  of  a 
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fixed-fixed  beam  and  a portal  frame,  and  Jensen  and  Sepulveda56  (1997)  minimized  the  weight 
of  a 25  bar  transmission  tower.  Fuzzy  set  theory  has  also  been  used  in  multidisciplinary 
optimization  by  Rao57  (1993),  who  designed  the  main  rotor  of  a helicopter  as  well  as  by  Wu 
and  Young58  (1996)  who  optimized  the  machine  room  layout  of  a ship.  Additionally,  Shih 
and  Chang59  (1995)  applied  multi-criteria  optimization  to  various  truss  examples  by 
considering  both  weight  and  displacement  as  objectives. 

Unfortunately,  designing  for  uncertainty  is  computationally  intensive,  typically 
requiring  an  order  of  magnitude  more  function  evaluations  than  a corresponding  deterministic 
design  problem.  In  reliability  optimization  using  probabilistic  methods,  approximation 
concepts  have  been  widely  used  to  reduce  the  computational  cost  e.g.,  Torng  and  Thacker60 
(1993),  Wang,  Grandhi  and  Hopkins26  (1994),  Sepulveda  and  Jensen61  (1995),  Wang  and 
Grandhi  28, 62  (1996),  Grandhi  and  Wang63  (1996)  all  used  various  approximation  concepts, 
although  none  of  them  considered  response  surface  approximations.  However,  Thacker  and 
Wu37  (1992)  used  response  surface  approximations  in  a probabilistic-based  reliability  analysis 
of  a turbine  blade  subject  to  creep  rupture,  while  Fox64  (1996)  investigated  various  aspects 
associated  with  using  response  surface  approximations  with  accurate  probabilistic  design. 

Compared  to  probabilistic  methods,  the  use  of  fuzzy  set  theory  to  model  uncertainty 
in  structural  optimization  problems,  is  a more  recent  development.  Sepulveda  and  Jensen65 
(1996)  and  Jensen  and  Sepulveda56  (1997)  reduced  the  computational  cost  associated  with 
performing  fuzzy  set  based  optimization  by  using  linear  Taylor  series  expansions  to  form  local 
approximations  that  are  applied  in  a sequential  manner.  In  fuzzy  set  theory,  the  high 
computational  cost  is  associated  with  the  evaluation  of  the  possibility  of  failure,  which 
requires  multiple  evaluations  of  a global  optimization  problem  that  is  defined  in  terms  of  the 


16 


uncertainty  problem  parameters.  Additionally,  the  possibility  of  failure  has  to  be  evaluated 
many  times  for  each  design  cycle.  In  terms  of  fuzzy  set  based  optimization,  response  surface 
approximations  have  the  advantage  that  the  number  of  function  evaluations  required  for 
constructing  response  surface  approximations  of  the  required  analysis  procedures,  are  typically 
less  than  the  number  of  function  evaluations  required  to  perform  a single  design  cycle. 
Furthermore,  once  the  response  surface  approximations  are  constructed  they  may  be  used 
repeatedly  without  requiring  any  new  function  evaluations,  thus  allowing  multiple 
optimizations  at  minimal  additional  computational  cost. 

In  the  present  dissertation,  the  possibility  of  failure  is  calculated  numerically  by 
making  use  of  the  vertex  method  (Dong  and  Shah66  1987),  which  is  based  on  the  definition  of  a 
fuzzy  function.  The  cost  of  evaluating  the  possibility  of  failure  is  reduced  by  replacing  costly 
finite  element  analyses  with  previously  developed  response  surface  approximations.  A 
response  surface  approximation  of  the  possibility  of  failure  is  constructed  from  numerical 
simulations,  which  involves  the  calculation  of  the  possibility  of  failure  for  different 
combinations  of  the  design  variables.  The  possibility  of  failure  response  surface 
approximation  is  thus  constructed  based  on  results  obtained  from  the  response  surface 
approximations  constructed  from  the  finite  element  data.  The  possibility  of  failure  response 
surface  approximation  simplifies  the  integration  of  the  analysis  code  and  the  optimization 
algorithm  and  filters  out  numerical  noise  in  the  approximate  response  function,  thus  enabling 
the  use  of  a derivative-based  optimization  algorithm. 

1 .6.  Objectives  of  Dissertation 

The  first  objective  of  the  present  dissertation  is  to  investigate  the  use  of  response 
surface  approximations  in  expensive  structural  optimization  problems,  and  to  develop  an 
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efficient  methodology  for  constructing  highly  accurate  response  surface  approximations  for 
use  in  structural  optimization  applications.  Even  though  the  use  of  response  surface 
approximations  in  design  optimization  has  several  advantages,  the  construction  of  response 
surface  approximations  is  not  a trivial  process  and  may  be  computationally  expensive  in  itself. 
Furthermore,  response  surface  approximations  are  a fairly  new  technology  in  the  area  of 
structural  optimization,  and  a need  exists  for  an  efficient  methodology,  tailored  specifically  to 
the  field  of  structural  optimization.  The  methodology  discussed  in  the  present  dissertation 
consists  of  a combination  of  both  existing  techniques  not  commonly  used  in  response  surface 
methodology  (e.g.,  dimensional  analysis  and  higher  order  polynomials  as  response  surface 
approximations)  as  well  as  some  existing  response  surface  methodology  techniques  not  often 
used  in  structural  optimization  applications  (e.g.,  stepwise  regression  and  detailed  error 
analysis). 

The  second  objective  of  this  dissertation  is  to  illustrate  the  advantages  of  using 
response  surface  approximations  in  expensive  structural  optimization  problems,  by 
considering  designing  for  uncertainty.  A stepped  plate  design  problem  is  considered  as  an 
illustrative  example,  and  response  surface  approximations  are  constructed  for  evaluating 
different  failure  criteria  of  the  plate  using  detailed  numerical  simulations  conducted  with  finite 
element  analyses.  Both  an  isotropic  as  well  as  a composite  laminated  plate,  where  the  change 
in  thickness  is  a result  of  an  internal  ply  drop  off,  are  considered.  The  problem  parameters  are 
uncertain  and  the  uncertainty  is  modeled  by  means  of  fuzzy  set  theory.  Response  surface 
approximations  are  used:  (1)  to  reduce  the  computational  cost  of  designing  for  uncertainty, 
(2)  to  integrate  the  analysis  code  and  the  optimization  algorithm,  and  (3)  to  filter  out 
numerical  noise  inherent  to  the  response  function. 
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1.7.  Outline  of  Dissertation 

Response  surface  methodology,  which  is  a key  aspect  of  the  present  dissertation,  is 
reviewed  in  Chapter  2.  The  review  concentrates  on  aspects  of  response  surface  methodology 
that  are  relevant  to  the  present  dissertation  and  the  use  of  response  surface  approximations  in 
structural  optimization  applications  in  general. 

Chapter  3 introduces  fuzzy  sets,  fuzzy  numbers  and  fuzzy  functions  and  describes  the 
relationship  between  the  membership  function  of  a fuzzy  number  and  the  possibility 
distribution  function.  Additionally,  the  vertex  method  is  introduced  as  a numerical  means  of 
calculating  membership  functions  for  fuzzy  functions  (a  function  of  fuzzy  numbers). 

In  Chapter  4 an  isotropic  stepped  plate  example  problem  and  its  associated  response 
surface  approximations  are  introduced  and  discussed  in  detail,  while  Chapter  5 summarizes  the 
optimization  results  obtained.  Similarly,  in  Chapter  6 a dropped  ply  composite  laminated 
plate  and  its  associated  response  surface  approximations  are  introduced  in  discussed  in  detail, 
with  the  corresponding  optimization  results  summarized  in  Chapter  7. 

Finally,  concluding  remarks  are  outlined  in  Chapter  8. 


CHAPTER  2 

RESPONSE  SURFACE  METHODOLOGY 

Response  surface  methodology  may  be  defined  as  “a  collection  of  statistical  and 
mathematical  techniques  useful  for  developing,  improving  and  optimizing  processes”  (Myers 
and  Montgomery6  1995,  p.  1)  or  as  “a  group  of  statistical  techniques  for  empirical  model 
building  and  model  exploration”  (Box  and  Draper7  1987,  p.  1).  The  methodology  was 
originally  developed  for  constructing  empirical  response  functions,  known  as  response  surface 
approximations  or  response  models,  based  on  a set  of  physical  experiments.  In  the  present 
dissertation,  the  physical  experiments  are  replaced  by  numerical  experiments  in  the  form  of 
finite  element  analyses.  Response  surface  methodology  typically  involves  the  identification  of 
a set  of  experimental  data  points  at  which  to  evaluate  the  response  function,  the  construction 
of  an  approximate  response  function  based  on  the  experimentally  obtained  data,  and  the 
evaluation  of  the  predictive  capabilities  of  the  approximate  response  function.  Finally,  the 
resulting  response  surface  approximations  are  used  in  optimization. 

2.1.  Response  Surface  Approximation  Construction 
A response  surface  approximations  is  an  approximate  relationship  between  a response 
1 1 (the  dependent  variable)  and  a vector  z of  k predictor  variables  (the  independent  variables). 
The  response  is  generally  obtained  from  experiments  (possibly  numerical  in  nature)  in  which 
case  7 denotes  the  mean  or  expected  response  value.  The  experimentally  obtained  response  y 
differs  from  the  expected  value  7 due  to  random  experimental  error  5.  The  relation  between  y 
and  the  predictor  variables  may  be  written  as 
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y = rj  + S = F(zl,z2,...,zk)  + S (2.1) 

where  F is  the  exact  dependence  of  77  on  z.  If  a response  surface  approximation  / is  used  to 
approximate  F,  Eq.  (2.1)  may  be  written  as 

y = f(zl,z2,...,zk)  + e (2.2) 

where  £ denotes  the  total  error,  which  is  the  difference  between  the  predicted  and  the 
expected  response  values,  and  includes  both  the  random  (variance)  error  8 and  the  modeling 
(bias)  error  (F-f).  The  response  surface  approximating  / is  generally  chosen  to  be  a low 
order  polynomial,  and  the  predicted  response  vector  y may  then  be  written  in  matrix  form  as 
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where  the  caret  symbol  denotes  predicted  values.  In  Eq.  (2.3),  n denotes  the  number  of  data 
points  considered,  and  p denotes  the  number  of  parameters  in  the  response  surface 
approximation.  Furthermore,  X is  a matrix  of  response  surface  approximation  parameters 
(usually  monomials)  evaluated  at  the  data  points,  and  the  vector  b contains  estimates  of  the 
unknown  coefficients  /?,  of  the  response  surface  approximation  f 

The  unknown  coefficients  ft,  of  the  response  surface  approximation  are  generally 
estimated  using  the  least  squares  procedure  (e.g.,  Burden  and  Faires2  1989,  pp.  427-439).  The 
least  squares  procedure  calculates  estimates  b,  of  the  actual  coefficients  ft,  in  such  a way  as  to 
minimize  the  sum  L of  the  squares  of  the  error  terms  as  follows: 

dL 


dp. 


= -2XTy  + 2XTXb  = 0 where  L = £f(2  = £(y, -/(z„z2  ..., zk))2  (2.4) 
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The  estimates  of  the  regression  coefficients  are  solved  from  Eq.  (2.4),  which  may  be  rewritten 
as  follows: 

b = (xTx)  *XTy  (2.5) 

The  difference  between  the  predicted  response  values  obtained  from  the  response 
surface  approximation  and  the  response  values  obtained  from  the  experiments  is  known  as  the 
residual  error,  defined  for  the  ith  data  point  as: 

ei  - y,  ~ yt  (2-6) 

2.2.  Statistical  Design  of  Experiments 

Statistical  design  of  experiments  is  an  integral  part  of  response  surface  methodology 
and  enables  the  designer  to  construct  response  surface  approximations  more  efficiently  by 
seeking  a small  set  of  optimal  data  points  at  which  to  evaluate  the  response  function.  This 
optimal  set  of  data  points  is  chosen  to  maximize  the  accuracy  by  minimizing  either  the 
variance  or  the  bias  error  associated  with  the  resulting  response  surface  approximation. 

In  cases  where  the  variance  error  is  dominant,  for  example  when  the  assumed  response 
surface  approximation  is  sufficient  to  accurately  approximate  the  actual  response  function  but 
the  data  are  subject  to  experimental  error,  it  is  appropriate  to  find  a set  of  optimal  data  points 
by  minimizing  the  variance  error.  For  regularly  shaped  design  spaces  there  are  a number  of 
standard  designs  (i.e.,  sets  of  data  points)  available  that  minimize  the  variance  error,  for 
example  the  Central  Composite  (Box  and  Wilson8  1951)  and  the  Box-Behnken  (Box  and 
Behnken9  1960)  designs.  For  an  irregularly  shaped  design  space,  the  designer  must  resort  to 
one  of  several  possible  computer  generated  designs  (e.g.,  Myers  and  Montgomery6  1995, 
pp.  364-366  and  Box  and  Draper7  1987,  pp.  489-501),  of  which  the  best  known  and  most 
widely  used  is  the  2>Optimality  design  criterion.  The  .D-Optimality  design  criterion 
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minimizes  a measure  of  the  variance  associated  with  the  estimates  of  the  coefficients  of  the 
response  surface  approximation,  by  maximizing  the  determinant  of  the  XrX  matrix.  When 
assuming  that  the  errors  e,  of  the  response  surface  approximation  are  independent  and  have  a 
normal  distribution  with  zero  mean  and  constant  variance,  the  determinant  of  the  XrX  matrix 
is  inversely  proportional  to  the  square  of  the  volume  of  the  confidence  region  associated  with 
the  regression  coefficients  b.  The  volume  of  the  confidence  region  provides  and  indication  of 
how  well  the  coefficients  of  the  response  surface  approximation  are  estimated,  with  a smaller 
volume  corresponding  to  better  estimation.  The  /9-Optimality  criterion  thus  selects  a set  of 
data  points  from  a candidate  set  of  data  points  to  maximize  the  determinant  of  the  XrX 
matrix,  and  the  best  /9-Optimal  design  is  obtained  from 

-f{^)  M 

where  £ represents  all  designs  in  the  set  of  possible  designs.  It  is  important  to  note  that 
variance-based  design  criteria  assume  that  all  of  the  error  associated  with  the  response  surface 
approximation  is  random,  and  do  not  address  bias  or  modeling  errors  resulting  from  an 
insufficient  response  surface  approximation. 

In  cases  where  the  bias  error  is  dominant,  for  example  when  low  order  polynomials 
are  used  to  approximate  a response  quantity  with  highly  non-linear  behavior,  it  may  be  more 
appropriate  to  find  a set  of  optimal  data  points  by  minimizing  the  bias  error  rather  than  the 
variance  error.  For  a regularly  shaped  design  space  and  a fixed  number  of  data  points,  closed 
form  solutions  exist  for  finding  the  minimum-bias  data  points  (e.g.,  Myers  and  Montgomery6 
1995,  pp.  406-421  and  Box  and  Draper7  1987,  pp.  437-442).  However,  for  an  irregularly  shaped 
design  space  no  closed  form  solution  exists  and  the  minimum-bias  data  points  must  be 
obtained  numerically  (e.g.,  Venter  and  Haftka11  1997). 
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Even  though  minimum-bias  statistical  designed  experiments  may  be  more  appropriate 
for  some  response  surface  approximations,  virtually  all  applications  within  the  engineering 
optimization  literature  make  successful  use  of  some  variance-based  statistical  design  of 
experiments.  In  the  present  dissertation  where  relatively  high  order  polynomials  (cubic  and 
quartic)  are  used  as  response  surface  approximations,  it  was  found  that  the  predictive 
capabilities  of  the  response  surface  approximations  were  not  significantly  improved  by  using  a 
minimum-bias  design  criterion.  Additionally,  it  is  difficult  to  implement  a minimum-bias 
design  for  an  irregular  design  space  with  a large  number  of  predictor  variables,  and  it  was  thus 
decided  to  use  the  /^-Optimality  criterion  whenever  statistical  design  of  experiments  was 
required.  In  all  cases,  the  25-Optimality  algorithm  provided  with  the  JMP67  (1995)  software 
package  was  used. 

2.3.  Elimination  of  Redundant  Parameters 
Since  the  exact  form  of  the  actual  response  function  is  rarely  known,  polynomial 
functions  are  generally  assumed  as  response  surface  approximations.  Due  to  the  general 
nature  of  such  a polynomial  response  function,  the  response  surface  approximation  will 
usually  include  redundant  parameters,  or  parameters  poorly  characterized  by  the  experiments. 
These  parameters  may  increase  the  error  of  the  response  surface  approximation,  thus 
decreasing  its  predictive  capabilities  (Myers  and  Montgomery6  1995,  p.  641,  and  Kaufman 
et  al.39  1996).  Eliminating  these  parameters  from  the  response  surface  approximation  is  an 
important  step  in  constructing  highly  accurate  response  surface  approximations. 

Several  methods  exist  for  selecting  the  best  subset  of  parameters  from  the  general 
polynomial  response  function  to  be  used  as  a final  response  surface  approximation.  Examples 
of  the  most  popular  procedures  are  the  so  called  “all-possible-regressions”  and  the  “stepwise 
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regression”  procedures  (e.g.,  Myers  and  Montgomery6  1995,  pp.  642-655  and  Ott68  1993, 
pp.  648-659).  The  all-possible-regressions  procedure  requires  that  all  possible  combinations  of 
parameters  from  the  general  response  surface  approximation  be  examined.  The  best  subset  of 
parameters  is  then  chosen  based  on  a suitable  criterion.  The  procedure  requires  the 
examination  of  2P  partial  response  surface  approximations  for  p parameters  in  the  general 
response  surface  approximation. 

The  stepwise  regression  procedure  was  developed  to  provide  an  alternative  to  the  high 
computational  cost  associated  with  the  all-possible-regressions  procedure,  especially  for 
response  surface  approximations  with  a large  number  of  parameters.  The  stepwise  regression 
procedure  represents  a family  of  procedures  that  examine  only  a small  number  of  partial 
response  surface  approximations  by  either  adding  or  deleting  one  parameter  at  a time  and  may 
be  classified  as  either  a forward,  a backward,  or  a mixed  stepwise  regression  procedure.  The 
forward  stepwise  regression  procedure  assumes  a response  surface  approximation  with  only 
one  parameter,  the  intercept.  Parameters  are  then  added  one  at  a time  by  selecting  the 
parameter  that  results  in  the  largest  increase  in  predictive  capabilities  of  the  current  response 
surface  approximation.  This  process  of  adding  parameters  is  repeated  until  a specified 
termination  criterion  is  satisfied.  In  contrast,  the  backward  stepwise  regression  procedure 
starts  with  a response  surface  approximation  that  includes  all  the  parameters  of  the  general 
response  surface  approximation.  From  this  response  surface  approximation,  one  parameter  is 
deleted  at  a time,  selecting  the  parameter  with  the  minimum  influence  on  the  predictive 
capabilities  of  the  resulting  response  surface  approximation.  Finally,  the  mixed  stepwise 
regression  procedure  consists  of  a combination  of  both  the  forward  and  backward  stepwise 
regression  procedures. 
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It  is  important  to  note  that  stepwise  regression  procedures  do  not  guarantee  the  best 
subset  of  parameters  and  must  be  used  with  caution.  However,  good  results  were  obtained  by 
using  the  mixed  stepwise  regression  procedure  provided  with  the  JMP67  (1995)  software 
package.  This  procedure  starts  with  an  initial  response  surface  approximation  that  includes  all 
the  parameters  of  the  general  response  surface  approximation  and  then  applies  backward 
stepwise  regression.  At  each  step  of  the  backward  stepwise  regression,  the  least  significant 
parameter  is  removed  from  the  current  response  surface  approximation.  Additionally,  at  each 
step  all  removed  parameters  are  re-examined  for  possible  inclusion  in  the  response  surface 
approximation.  The  procedure  may  be  summarized  as: 

1.  Start  with  the  full  response  surface  approximation  as  the  initial  response  surface 
approximation. 

2.  Apply  backward  stepwise  regression,  thus  deleting  the  least  significant  parameter  from 
the  current  response  surface  approximation. 

3.  Apply  forward  stepwise  regression  to  all  of  the  deleted  parameters.  This  procedure 
evaluates  if  any  of  the  deleted  parameters  has  a significant  influence  on  the  predictive 
capability  of  the  current  response  surface  approximation.  The  criterion  used  to 
determine  when  to  include  a parameter  is  generally  more  stringent  than  the  criterion 
used  to  remove  a parameter. 

4.  Go  to  step  2 and  repeat  until  some  termination  criterion  is  satisfied. 

In  the  present  dissertation,  Mallow’s  Cp  statistic  was  used  as  the  termination  criterion 
to  identify  the  best  reduced  response  surface  approximation  from  the  subset  of  reduced 
response  surface  approximations  provided  by  the  stepwise  regression  procedure.  The  Cp 
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statistic  addresses  the  prediction  capabilities  of  the  response  surface  approximation  by 
estimating  a total  error.  This  total  error  consists  of  two  parts,  one  resulting  from  modeling 
error  (bias)  and  the  other  from  noise  (variance  error).  An  under  specified  response  surface 
approximation  (not  enough  parameters)  increases  that  part  of  the  total  error  resulting  from 
bias,  while  an  over  specified  response  surface  approximation  (too  many  parameters)  increases 
that  part  of  the  total  error  resulting  from  variance  (Ott68  1993,  p.  653).  The  Cp  statistic  is 
defined  as 
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where  SSEP  is  the  sum  of  the  squares  of  the  n errors  for  a response  surface  approximation  with 
p parameters  (including  the  intercept  or  constant  parameter)  and  MSE  is  the  mean  sum  of 
squares  of  the  error  terms  obtained  from  the  response  surface  approximation  with  all 
parameters  included.  Based  on  the  Cp  statistic,  the  best  response  surface  approximation  in  a 
set  of  candidate  response  surface  approximations  is  the  response  surface  approximation  with 
the  lowest  Cp  value,  and  a good  response  surface  approximation  should  at  least  have  a Cp 
value  close  to  p. 


2.4.  Evaluation  of  Predictive  Capabilities 
Optimization  has  the  general  tendency  of  exploiting  weaknesses  in  the  formulation  of 
the  response  function.  In  terms  of  response  surface  approximations,  this  means  that  the 
optimizer  would  tend  to  drive  the  optimum  design  to  a region  where  the  approximations  are 
inaccurate.  Highly  accurate  response  surface  approximations  are  thus  a requirement  for 
structural  optimization  applications,  and  it  is  important  to  evaluate  the  predictive  capabilities 
of  the  approximations.  In  the  present  dissertation,  a detailed  error  analysis  using  two  data  sets 
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and  various  statistical  criteria  were  used  to  evaluate  the  predictive  capabilities  of  the  response 
surface  approximations.  Data  Set  1 contains  all  data  points  used  in  constructing  the  response 
surface  approximations,  while  Data  Set  2 is  an  independent  data  set. 

For  Data  Set  1,  the  coefficient  of  determination  (R2)  statistic,  the  adjusted  R2  ( Adj-R I2) 
statistic,  the  percent  root  mean  square  error  ( %RMSE ),  and  the  percent  root  mean  square  error 
based  on  the  predicted  sum  of  squares  ( PRESS)  statistic  ( %RMSEpress ) are  calculated  (Myers 
and  Montgomery6  1995,  pp.  28-47).  The  coefficient  of  determination  (R2  ) statistic  is  the 
proportion  of  the  variability  in  the  response  data  that  is  accounted  for  by  the  response  surface 
approximation  and  is  defined  as 


R2  =1 
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The  coefficient  of  determination  has  a value  between  0 and  1.  Unfortunately,  a value  close 
to  1 does  not  necessarily  imply  a good  response  surface  approximation  because  additional 
parameters  added  to  a response  surface  approximation  will  always  increase  the  R2  value 
without  necessarily  increasing  the  predictive  capabilities  of  the  response  surface 
approximation.  The  adjusted  coefficient  of  determination  (Adj-R? ) statistic  is  an  alternative 
measure  of  the  explained  variability  that  is  often  used  and  has  the  desirable  property  that  its 
value  does  not  necessarily  increase  when  adding  additional  parameters  to  a response  surface 
approximation.  When  a large  difference  exists  between  the  R2  and  Adj-R2  values,  there  exists  a 
good  chance  that  the  response  surface  approximation  contains  insignificant  parameters.  The 
Adj-R2  statistic  is  defined  as: 
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(2.11) 


The  percentage  errors  are  defined  in  relation  to  the  average  absolute  value  y of  the 
measured  response  data,  where 


y=-tM 


(2.12) 


and  the  %RMSE  is  defined  as 


o%RMSE  = ~\-±(yi-y1)2  (2.13) 

y 

where  % = (n  - p)  for  Data  Set  1,  and  % = n for  Data  Set  2.  For  Data  Set  1,  the  ( n - p)  factor 
is  used  to  obtain  an  unbiased  estimator  of  the  mean  square  error  of  the  response  surface 
approximation,  thus  compensating  for  the  fact  that  the  error  is  calculated  at  the  data  points 
used  to  construct  the  response  surface  approximation.  Even  though  Eq.  (2.13)  is  based  on  an 
unbiased  estimator  of  the  mean  square  error,  the  %RMSE  value  may  be  overly  optimistic  in 
estimating  the  predictive  capabilities  of  the  response  surface  approximation.  A more  realistic 
estimate  of  the  predictive  capabilities  of  a response  surface  approximation  is  obtained  when 
defining  the  %RMSE  based  on  the  PRESS  statistic.  The  PRESS  statistic  is  calculated  by 
selecting  a data  point,  say  data  point  /.  The  response  surface  approximation  obtained  from 
the  remaining  (n  - 1)  data  points  is  then  used  to  predict  the  response  at  the  withheld  data 
point,  denoted  by  y(i).  The  prediction  error  at  the  withheld  data  point  e<t)  may  then  be 

written  as 

e(l)  =y<-y(l)  (2-14) 

and  is  referred  to  as  the  Ith  PRESS  residual.  This  procedure  is  then  repeated  for  all  data  points 
and  the  resulting  PRESS  residuals  are  summed  to  obtain  the  PRESS  statistic  as  follows: 
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PRESS  = £ef0  =£(?,-  yj  (2.15) 

i=\  i=l 

However,  from  (Myers  and  Montgomery6  1995,  p.  43-47)  the  /*  PRESS  residual  may  also  be 


written  as 


e(.t)  ~ 


1-/7,. 


(2.16) 


where  /?,,  denotes  the  diagonal  elements  of  the  matrix  H,  which  is  define  by: 

H = X(XrX)“‘Xr  (2.17) 

The  formulation  of  shown  in  Eq.  (2.16)  is  significant  since  it  allows  the  calculation  of  the 
PRESS  statistic  from  a single  linear  regression,  by  using  the  relationship: 


PRESS  = £ 


( V 
e, 


(2.18) 


The  %RMSE  based  on  the  PRESS  statistic  ( %RMSEPREss ) is  then  defined  as: 


% RMSEP 


y V n 


100  - PRESS 


(2.19) 


Apart  from  the  %RMSE,  the  percent  average  error  ( %AvgErr  ) and  the  percent  maximum 
error  ( %MaxErr ) are  also  calculated  for  Data  Set  2 as  follows: 


% AvgErr  = -y\ 

ny  ti 


VoMaxErr  = Max  \y  - y. 

' y 


(2.20) 


(2.21) 


2.5.  Outliers.  Leverage  and  Influential  Data  Points 
Apart  from  evaluating  the  predictive  capabilities  of  the  constructed  response  surface 
approximations,  it  is  also  important  to  evaluate  the  quality  of  the  data  points  used  to  construct 
the  response  surface  approximations.  An  evaluation  of  the  data  points  is  required  since  it  is 
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possible  for  a small  set  of  data  points,  known  as  influential  data  points,  to  have  a 
disproportionately  high  influence  on  the  resulting  response  surface  approximation.  An 
influential  data  point  is  defined  as  a data  point  which,  if  removed,  the  resulting  response 
surface  approximation  changes  considerably.  Influential  data  points  are  typically  data  points 
with  high  leverage  and  an  outlier  tendency. 

The  leverage  of  a data  point  is  related  to  the  location  of  that  data  point  relative  to  the 
centroid  of  all  data  points  used  to  construct  the  response  surface  approximation.  Data  points 
located  far  from  the  centroid  have  high  leverage.  The  diagonal  elements  htt  of  the  matrix  H 
Eq.  (2.17)  are  useful  in  determining  the  leverage  of  data  points.  The  diagonal  elements  /?,,  are 
bounded  by  0 <hu<  1 with  an  average  value  of  p/n.  According  to  Myers  and  Montgomery6 
(1995,  p.  49)  a rough  guideline  is  that  a ha  value  greater  than  2p/n  would  indicate  that  the  Xth 
data  point  is  a high  leverage  data  point. 

An  outlier  is  a data  point  with  a large  residual  error  e,  and  may  be  a result  of  either 
inaccurate  experimental  data  or  an  inadequate  response  surface  approximation.  Outliers  are 
typically  detected  by  using  the  ^-Student  residual  statistic  t„  which  is  obtained  by  scaling  the 
residual  e*  of  the  ith  data  point  by  (Myers  and  Montgomery6  1995,  pp.  43-48) 


t.  = 


and 


4 = 


(n  - p)MSE  - e]  7(1  - hu ) 
n- p-l 


(2.22) 


where  % is  an  estimate  of  the  variance  of  the  error  terms  of  the  data  set  with  the  ith  data  point 
removed.  Compared  to  the  unsealed  residuals,  the  7?-Student  residuals  have  constant  variance 
Var(t)  = l regardless  of  the  position  of  the  i‘h  data  point  and  are  more  sensitive  to  influential 
data  points.  Any  value  outside  the  interval  -3  < ?,  < 3 should  be  considered  as  possible 


outliers. 
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Influential  data  points  are  typically  difficult  to  detect  since  their  residuals  may  be 
artificially  small  due  to  the  large  influence  these  data  points  have  on  the  resulting  response 
surface  approximation.  The  Cook's  D statistic  is  used  to  identify  influential  data  points  and  is 
defined  as  (Myers  and  Montgomery6  1995,  pp.  49-50): 

D (b(l)  ~b)rXrX(b(<)  -b) 

' pMSE 

The  Cook’s  D statistic  is  a measure  of  the  squared  distance  between  the  vector  of  least  squares 
estimated  coefficients  b obtained  from  all  of  the  n data  points  and  the  vector  of  least  squares 
estimated  coefficients  b(/)  obtained  by  deleting  the  ith  data  point  from  the  original  set  of  data 
points.  According  to  Myers  and  Montgomery6  (1995,  pp.  49-50),  data  points  with  a Cook’s  D 
statistic  larger  than  1 is  considered  as  influential  data  points.  Influential  data  points  may  also 
be  detected  by  using  robust  regression  procedures  which  are  designed  to  dampen  the  effect 
that  influential  data  points  may  have  on  a response  surface  approximation  obtained  from  the 
normal  least  squares  procedure.  Robust  regression  is  typically  applied  using  an  iteratively 
reweighted  least  squares  procedure  that  assigns  weights  to  all  data  points  used  to  construct  the 
response  surface  approximation.  Influential  data  points  are  assigned  a low  weight,  thus 
increasing  the  associated  residuals  (as  compared  to  those  obtained  from  the  normal  least 
squares  procedure)  and  making  it  much  easier  to  identify  influential  data  points.  A detailed 
discussion  of  robust  regression  techniques  are  beyond  the  scope  of  this  dissertation  and  the 
interested  reader  is  referred  to  Myers  and  Montgomery6  (1995,  pp.  667-673)  and  Box  and 
Draper7(1987,  pp.  85-90)  as  an  introduction  with  several  additional  references.  When 
influential  data  points  are  detected,  these  data  points  should  be  carefully  investigated  to 
determine  if  they  are  influential  as  a result  of  inaccurate  experimental  data  or  due  to  an 
inadequate  response  surface  approximation,  before  taking  appropriate  action. 
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In  order  to  evaluate  the  need  for  an  extensive  investigation  to  identify  influential  data 
points  in  the  data  sets  of  the  present  dissertation,  a few  representative  sample  data  sets  were 
investigated  first.  These  data  sets  were  investigated  using  the  /?-Student  and  Cook’s  D 
statistics  as  well  as  an  iteratively  reweighted  least  squares  procedure  to  ensure  that  the  data 
points  were  of  good  quality.  It  was  found  that  since  accurate  finite  element  analyses  and 
appropriate  high  order  polynomials  were  used  in  constructing  the  response  surface 
approximations,  no  influential  data  points  were  present.  Note  that  the  use  of  a second, 
independent  data  set  in  calculating  the  predictive  capabilities  of  the  response  surface 
approximations  (Section  2.4)  also  protects  against  possible  influential  data  points  and  was  the 
only  diagnostic  considered  to  ensure  data  points  of  good  quality  throughout  the  rest  of  this 
dissertation. 


2.6.  Optimization 

One  of  the  most  attractive  features  of  using  response  surface  approximations  in 
optimization  is  that  due  to  the  least  squares  procedure  and  the  use  of  polynomial  functions  as 
response  surface  approximations,  the  resulting  approximate  response  function  is  a smooth 
function.  Having  a smooth  response  function  is  important  since  most  optimization 
algorithms  have  problems  dealing  with  non-smooth  response  functions  because  the 
optimization  algorithm  easily  converges  on  a local  minimum.  In  the  case  of  numerical 
experiments,  the  non-smooth  functional  behavior  is  a result  of  numerical  noise,  which  creates 
spurious  local  minima  that  easily  trap  the  optimization  algorithm  and  make  it  difficult  to  find 
the  best  designs  (e.g.,  Giunta3  et  al.  1994). 

When  replacing  the  noisy  response  function  obtained  from  the  numerical  experiments 
with  a smooth  response  surface  approximation,  one  effectively  filters  out  the  numerical  noise, 
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thus  allowing  the  use  of  derivative-based  optimization  algorithms.  Derivative-based 
optimization  algorithms  are  efficient  algorithms  and  are  readily  available  commercially.  In 
the  present  dissertation,  both  the  generalized  reduced  gradient  algorithm  provided  with 
Microsoft  Excel  Version  8.069  (1997)  and  the  sequential  linear  programming  (SLP)  procedure 
provided  with  the  Design  Optimization  Tools70  (DOT)  (1995)  software,  were  used  to  solve  the 
optimization  problems.  Microsoft  Excel  Version  8.069  (1997)  provided  an  easy  to  use, 
interactive  environment.  However,  since  each  optimization  problem  is  solved  interactively, 
this  environment  is  not  suitable  for  performing  a large  number  of  optimizations.  In  contrast, 
the  Design  Optimization  Tools70  (1995)  software  is  supplied  as  a number  of  FORTRAN71 
(1995)  subroutines.  These  subroutines  can  easily  be  called  from  a simple  FORTRAN71  (1995) 
program  to  perform  a large  number  of  optimizations  with  no  further  interaction  required 
from  the  designer. 


CHAPTER  3 
FUZZY  SET  THEORY 


Fuzzy  set  theory  was  introduced  by  Zadeh53  (1965)  as  a mathematical  tool  for  the 
quantitative  modeling  of  uncertainty,  and  makes  use  of  fuzzy  numbers  to  represent  uncertain 
problem  parameters.  In  order  to  represent  an  uncertain  parameter  as  a fuzzy  number,  the 
designer  only  needs  to  specify  the  range  of  uncertainty  and  a membership  function  that 
denotes  the  possibility  of  occurrence  of  an  element  in  the  specified  range.  Membership 
functions  are  generally  constructed  subjectively,  based  on  expert  opinion.  Fuzzy  set  theory 
caters  to  worst  case  scenarios  and  thus  compensates  for  the  fact  that  the  uncertainty  is 
modeled  based  on  subjective  opinions  and  assumptions  (Maglaras  et  al.52  1997). 

3.1.  Fuzzy  Sets,  Fuzzy  Numbers  and  Fuzzy  Functions 
Fuzzy  set  theory  presents  a methodology  for  the  mathematical  modeling  of 
uncertainty.  In  contrast  to  classical  (or  crisp)  set  theory  where  an  element  either  completely 
belongs  to  a specific  set  or  not  at  all,  fuzzy  set  theory  makes  use  of  membership  functions  to 
denote  the  degree  to  which  an  element  belongs  to  a fuzzy  set. 

Any  crisp  set  has  a characteristic  equation  that  assigns  a value  of  either  0 or  1 to  each 
member  of  the  universal  set,  where  0 is  generally  used  to  denote  non-membership  and  1 to 
denote  membership.  An  example  is  a set  B defined  such  that  x is  any  real  number  between  2 
and  3,  which  may  be  written  as: 

5 = {x|2<x<3,xe9l}  (3.1) 

The  characteristic  equation  for  set  A is  then  defined  as: 
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J 1 for  x e B 
^ B [0  for  x g B 

Thus,  for  the  crisp  set  B the  characteristic  equation  maps  elements  of  the  set  91  to  elements  of 
the  set  {0,1}  as  follows: 

^:5R->{0,1}  (3-3) 

A fuzzy  set  is  defined  by  generalizing  the  characteristic  equation  to  assign  a grade  of 
membership,  ranging  between  0 and  1,  to  each  element  of  the  universal  set.  In  fuzzy  set 
theory,  this  generalized  characteristic  equation  is  referred  to  as  a membership  function.  The 
notation  used  to  denote  a membership  function  is  as  follows: 

M(x):Y  ->  [0,1]  (3.4) 

Where  M denotes  the  membership  function,  mapping  the  elements  of  the  universal  set  X to 
the  real  interval  [0,  1].  In  the  present  dissertation  the  same  symbol,  a bold  face  capital  letter,  is 
used  to  denote  both  the  fuzzy  set  and  its  membership  function.  Since  each  fuzzy  set  is 
completely  and  uniquely  defined  by  only  one  particular  membership  function,  no  ambiguity 
results  from  the  double  use  of  the  same  symbol. 

Fuzzy  sets  may  be  represented  numerically  by  making  use  of  a level  cuts.  An  a level 
cut  is  defined  as  the  real  interval  where  the  membership  function  is  larger  than  a given  value  a 
(Klir  and  Yuan72  1995,  p.  19)  and  may  be  written  mathematically  for  a generic  fuzzy  set  B as 
follows: 

“A  = {r|B(x)>a}  (3.5) 

Equation  (3.5)  is  shown  graphically  in  Fig.  (3.1),  where  it  is  assumed  that  B has  a 
triangular  and  symmetric  membership  function.  Also  shown  in  Fig.  (3.1)  are  the  lower  and 
upper  bounds  of  the  a level  cut,  ab\  and  ab2. 


36 


Figure  3.1:  Triangular  and  symmetric  membership  function  with  a level  cut  shown. 

A fuzzy  number  is  defined  as  a fuzzy  set  that  is  both  normal  and  convex  (Klir  and 
Yuan72  1995,  p.  97).  A normal  fuzzy  set  is  one  for  which  the  maximum  membership  function 
is  equal  to  1,  while  a convex  fuzzy  set  is  one  for  which  all  possible  a level  cuts  are  convex. 
The  fuzzy  set  B shown  in  Fig.  (3.1)  is  thus  a fuzzy  number.  In  fact,  the  triangular  and 
symmetric  membership  function  is  most  often  used  to  represent  fuzzy  numbers,  mainly  due 
to  its  simplicity.  As  shown  in  Fig.  (3.1)  any  triangular  fuzzy  number  (not  necessarily 
symmetric)  may  be  represented  by  only  three  variables:  xL\xc\  and  xR. 

A fuzzy  function  Y is  a function  of  fuzzy  variables  X,  and  may  be  written  as 

Y = Y(X1,X2,...,X„)  (3.6) 

for  the  case  where  n fuzzy  variables  are  considered.  The  following  properties  of  fuzzy 
functions  are  summarized  and  proved  by  Klir  and  Yuan72  (1995,  pp.  105-109): 

1.  When  all  of  the  fuzzy  variables  of  a fuzzy  function  are  continuous  fuzzy  numbers,  the 
fuzzy  function  itself  is  also  a continuous  fuzzy  number. 
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2.  When  all  of  the  fuzzy  variables  of  a fuzzy  function  are  fuzzy  numbers,  the  a level  cut 
of  a fuzzy  function  aY  may  be  written  in  terms  of  the  a level  cuts  of  its  fuzzy  variables 
aX,  as  follows 


aY 


Y(aXvaX2,...,°Xn) 

min{Y (aXx ,aX2,.. .°Xn )} max(Y{°Xx  °X2,.. .“Xn )) 
. aR  aR 


(3.7) 


where  aR  denotes  the  ^-dimensional  rectangle,  formed  by  the  a level  cuts  of  the  n 
fuzzy  numbers. 


3.2.  The  Vertex  Method 

Based  on  the  properties  of  a fuzzy  function  (see  Eq.  3.7),  Dong  and  Shah66  (1987) 
introduced  the  vertex  method  for  evaluating  the  upper  and  lower  bounds  of  aY  when  all  the 
fuzzy  variables  of  Y are  fuzzy  numbers.  When  aY  has  no  global  extreme  points  in  the  interior 
(including  the  boundaries)  of  the  ^-dimensional  rectangle  formed  by  the  a level  cuts  of  the  n 
fuzzy  numbers,  the  vertex  method  yields 


aY  = 


mm 

J 


in(Y(acJ)]mqx(YCcJ)) 


V y = 1,2" 


(3.8) 


where  aCj  =(aXvaX2,...,aXn)J  denotes  the  coordinates  of  the  2”  vertices  of  the  M-dimensional 

rectangle  obtained  from  the  fuzzy  numbers  for  a specific  a value.  When  the  fuzzy  function 
has  global  extreme  points  in  this  ^-dimensional  rectangle  (including  the  boundaries),  the 
function  evaluations  at  these  extreme  points  must  be  included  when  evaluating  the  minimum 
and  maximum  operators  of  Eq.  (3.8).  Note  that  since  all  possible  combinations  of  the  a level 
cut  end  points  of  the  fuzzy  numbers  are  considered,  the  evaluation  of  Eq.  (3.8)  generally 
requires  a large  number  of  function  evaluations  and  is  computationally  intensive. 
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In  order  to  ensure  that  no  global  extreme  points  exist  within  the  design  space 
(including  the  boundaries),  except  at  the  vertices,  the  sequential  linear  programming  procedure 
provided  with  the  Design  Optimization  Tools  (DOT)70  (1995)  software  was  used.  A large 
number  of  derivative-based  optimizations  from  different  starting  points,  randomly  distributed 
throughout  the  design  space,  were  performed  for  a small  number  of  representative  a level  cut 
values.  The  results  obtained  from  these  optimizations  were  then  investigated  to  ensure  that 
no  global  extreme  points  existed  at  locations  other  than  the  vertices  of  the  design  space. 

3.3.  Possibility  Measure  and  Possibility  Distribution 
A fuzzy  number  may  also  be  considered  as  the  trace  of  a possibility  measure  n on  the 
singletons  (single  elements)  x of  the  universal  set  X (Dubois  and  Prade73  1988,  pp.  13-17). 
When  a possibility  measure  defined  on  the  unit  interval  is  considered,  its  possibility 
distribution  n can  then  be  interpreted  as  the  membership  function  of  a fuzzy  number  B 
describing  the  event  that  n focuses  on,  which  may  be  written  mathematically  as: 

n({x})  = ^(x)  = B(x)  V xeX  (3.9) 

The  possibility  of  an  event  is  calculated  from  the  possibility  distribution  function  and 
is  related  to  the  maximum  value  of  the  possibility  distribution  function.  The  calculation  of 
the  possibility  of  an  event  may  be  illustrated  with  a simple  example  problem  where  a real 
number  x that  is  approximately  equal  to  2 is  considered.  This  number  may  be  represented  by 
a fuzzy  set  B with  a symmetric  and  triangular  membership  function  (see  Fig.  3.2)  where  it  is 
assumed  that  real  numbers  between  1.5  and  2.5  have  non-zero  possibility  of  being 
approximately  equal  to  2. 

Recall  that  according  to  Eq.  (3.9)  the  membership  function  of  Fig.  (3.2)  is  equivalent  to 
the  possibility  distribution  for  the  event  that  a given  real  number  x is  approximately  equal 
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Figure  3.2:  Possibility  distribution  for  x approximately  equal  to  2. 

to  2.  The  possibility  that  the  value  of  x would  fall  in  a specified  range,  say  from  1.5  to  1.75,  is 
equal  to  the  maximum  value  of  the  possibility  distribution  function  (i.e.,  the  membership 
function)  over  that  range,  in  this  case  0.5.  Similarly,  the  possibility  that  the  value  of  x would 
fall  in  the  range  from  1.5  to  2.0  is  equal  1.0,  which  is  the  same  possibility  that  the  value  of  x 
would  fall  in  the  range  from  1.75  to  2.0.  Additionally,  the  possibility  that  the  value  of  x is 
exactly  equal  to  a specific  value,  say  1.75,  is  equal  to  the  value  of  the  possibility  distribution  at 
that  value,  in  this  case  equal  to  0.5.  Similarly,  the  possibility  that  the  value  of  x is  equal  to  2.0 
would  thus  be  equal  to  1 .0. 

Calculating  the  possibility  of  an  event  is  significantly  different  from  calculating  the 
probability  of  the  same  event.  The  probability  of  an  event  is  calculated  from  a probability 
distribution  function  and  is  related  to  the  area  under  the  probability  distribution  function.  To 
illustrate  calculating  the  probability  of  an  event,  an  uniform  probability  distribution  function 
may  be  used  to  represent  the  real  number  x that  is  approximately  equal  to  2.  The  uniform 
probability  distribution  function  indicates  that  the  value  of  x has  an  equal  probability  of  being 
equal  to  any  real  number  in  the  range  1.5  to  2.5,  and  is  shown  in  Fig.  (3.3). 
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Figure  3.3:  Uniform  probability  density  function  for  x approximately  equal  to  2. 

Note  that  a probability  distribution  function  is  always  constructed  so  that  the  total  area  under 
the  graph  is  equal  to  1.  The  probability  that  the  value  of  x would  fall  in  a specified  range,  say 
from  1.5  to  1.75,  would  then  be  equal  to  the  area  under  that  part  of  the  probability 
distribution  function,  in  this  case  0.25.  Similarly,  the  probability  that  the  value  of  x would  fall 
in  the  range  from  1.75  to  2.0  is  also  equal  to  0.25  and  in  the  range  from  1.5  to  2.0  is  equal 
to  0.5.  Additionally,  the  probability  that  the  value  of  x would  be  equal  to  a specific  number, 
for  example  1.75  or  2.0,  is  equal  to  zero  since  the  area  under  that  part  of  the  graph  (a  point)  is 
equal  to  zero. 

The  difference  between  the  probability  and  possibility  of  an  event  stems  from  the  fact 
that  probability  and  possibility  are  two  different  measures  of  uncertainty  (Dubois  and  Prade73 
1988,  pp.  1-4).  Probability  theory  is  concerned  with  the  frequency  of  the  occurrence  of  an 
event.  For  our  example  problem,  the  event  may  be  drawing  a card  from  a bag  with  an  infinite 
number  of  cards  each  denoting  a real  number  between  1.5  and  2.5.  Thus,  the  probability  of 
drawing  a card  that  has  a value  exactly  equal  to  a specific  real  number  (say  1.75)  would  be 
equal  to  0,  while  drawing  a card  with  a number  in  the  range  1.5  to  1.75  would  be 
equal  to  0.25. 
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In  contrast,  the  possibility  of  an  event  is  a measure  of  how  plausible  that  event  is,  in 
other  words  how  surprised  you  would  be  if  you  knew  x is  approximately  equal  to  2.0  and  the 
actual  value  turns  out  to  be  equal  to  1.75  for  example.  A more  detailed  discussion  on  the 
differences  between  probability  theory  and  possibility  theory  are  provided  in  Dubois  and 
Prade73  (1988,  pp.  1-13)  and  Klir  and  Yuan72  (1995,  pp.  200-208). 

Finally,  for  calculating  the  possibility  of  failure  it  is  often  required  to  compare  a crisp 
number  with  a fuzzy  number.  The  possibility  measure  of  a crisp  number  x being  smaller  or 
equal  to  a fuzzy  number  B is  defined  by  Dubois  and  Prade73  (1988,  pp.  99-101)  as  follows: 

nB([x,+oo))  = supB(y)  Vx  P-10) 


y^x 


The  possibility  distribution  function  /rB  corresponding  to  the  possibility  measure  of 
Eq.  (3.10)  is  shown  graphically  in  Fig.  (3.4)  for  the  general  case  where  B has  a non-linear 
membership  function. 


Figure  3.4:  Possibility  distribution  for  B > x with  non-linear  B(x). 

Using  the  definition  of  Eq.  (3.10),  the  possibility  distribution  of  failure  ;r(P_P/)  is 

obtained  from  the  fuzzy  function  (P  - P/)  which  contains  the  fuzzy  numbers  P (the  applied 
load)  and  P/  (the  failure  load  of  the  structure)  as  fuzzy  variables.  The  possibility  of  failure  is 


then  defined  as: 
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n(p-p,)([°,-HX>))  = sup(p  - p/Xy)  (3*n) 

J y^O 

In  the  present  dissertation,  the  programming  languages  FORTRAN71  (1995)  and 
JAVA74  (1996)  were  used  to  develop  software  for  calculating  the  possibility  of  failure  where  all 
uncertain  input  parameters  were  represented  by  triangular  and  symmetric  membership 
functions. 


CHAPTER  4 

ISOTROPIC  PLATE  RESPONSE  SURFACE  CONSTRUCTION 

Designing  for  uncertainly  of  an  isotropic  plate  with  a change  in  thickness  across  its 
width  is  the  first  example  problem  that  will  be  considered.  The  uncertainty  considered  is 
typical  of  the  aircraft  industry  where  structures  are  often  designed  that  will  be  built  well  into 
the  future  from  materials  not  yet  available.  Oftentimes  these  materials  are  not  yet  available 
and  the  designer  is  forced  to  deal  with  uncertain  material  properties.  In  addition,  the  cost  of 
manufacturing  is  also  uncertain.  However,  unlike  the  material  properties,  the  designer  has 
some  control  over  the  manufacturing  cost  which  is  closely  linked  to  the  required  tolerances  in 
geometry.  For  such  design  problems,  little  information  about  the  uncertainty  is  known  and 
the  uncertainty  is  typically  based  on  expert  opinion  and  assumptions  made  by  the  designer. 
Fuzzy  set  theory  is  ideally  suited  for  this  class  of  problems  and  is  capable  of  modeling 
uncertainty  using  limited  data  available  to  the  designer. 

4.1.  Problem  Description 

An  isotropic  plate  with  a change  in  thickness  in  the  form  of  a linear  ramp  (see  Fig.  4.1) 
is  the  first  design  problem  considered  in  the  present  dissertation.  The  plate  is  simply 
supported  on  two  of  its  edges,  free  on  the  other  two  edges,  and  subjected  to  a uniformly 
distributed  load,  applied  on  the  two  simply  supported  edges.  Three  non-dimensional 
parameters  A,  P and  y specify  the  geometry  and  location  of  the  change  in  thickness  (see 
Fig.  4.2).  These  variables  vary  in  the  following  intervals:  A e [0,1],  P e [-0.5, 0.5]  and  ye  [0,1]. 
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Figure  4.1: 


Three  dimensional  perspective. 


Figure  4.2:  Cross-section  model  with  non-dimensional  variables  shown. 

Yielding  is  considered  as  a failure  mechanism  in  the  design,  while  the  plate  is  also 
designed  to  be  buckling  resistant.  The  Von  Mises  failure  criterion  is  used  to  calculate  the 
failure  load  of  the  plate  due  to  yielding.  For  the  buckling  load  constraint,  the  design  fails  if 
the  applied  load  exceeds  the  buckling  load.  Response  surface  approximations  for  the  stress 
distribution  near  the  re-entrant  corner  and  for  the  buckling  load  were  constructed  in  terms  of 
the  geometry  and  location  of  the  change  in  thickness.  Numerical  experiments  were  conducted 
using  finite  element  models  of  the  plate  with  MSC/NASTRAN  Version  6875  (1994).  A cross- 
section  of  the  plate  was  used  to  model  the  stress  distribution  near  the  re-entrant  corner,  using 
four-node,  isoparametric,  plane  strain  elements.  All  models  had  a uniform  finite  element 
mesh,  typically  having  1,800  elements  in  the  x-direction  and  9 elements  in  the  z-direction.  The 
number  of  elements  varied  slightly  from  model  to  model.  A schematic  representation  of  the 
finite  element  model  used  is  shown  in  Fig.  (4.3). 
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Figure  4.3:  Schematic  finite  element  model  used  to  determine  the  stress  distribution  near 

the  re-entrant  corner  (actual  model  has  about  16,000  elements). 

For  the  buckling  response,  constant  thickness,  four-node,  isoparametric  plate  bending 
elements  were  used  to  construct  a two-dimensional  twenty  element  by  twenty  element  finite 
element  model  similar  in  shape  to  a plan  view  of  Fig.  (4.1),  as  shown  in  Fig  (4.4).  The 
thickness  of  the  elements  in  the  ramp  region  is  the  average  thickness  over  the  element.  The 
eccentricity  of  the  mid-plane  was  found  to  have  an  insignificant  effect  on  the  Euler  buckling 
load  value  and  was  ignored  in  the  analysis.  The  eccentricity  is  expected  to  cause  large 
deflections  when  the  applied  load  approaches  the  buckling  load  value,  which  usually  results  in 
stress  failure. 


Figure  4.4:  Schematic  finite  element  model  used  to  determine  the  buckling  load  (actual 

model  has  400  elements). 

4.2.  Dimensional  Analysis 


By  reducing  the  number  of  variables  in  a response  surface  approximation,  the  accuracy 
and  the  usefulness  of  the  response  surface  approximation  are  increased.  The  number  of 
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variables  in  the  response  surface  approximations  of  the  present  example  was  reduced  by  using 
variables  that  are  intrinsic  to  the  problem,  identified  by  performing  a dimensional  analysis  of 
the  governing  equations. 

4.2.1. Stress  Distribution  near  the  Re-entrant  Corner  of  the  Plate 

The  theoretical  value  of  the  stress  at  the  re-entrant  corner  is  infinite.  However,  due  to 
yielding  and  the  impossibility  of  manufacturing  an  absolutely  sharp  corner,  the  actual  stress  at 
the  corner  will  be  finite.  To  account  for  yielding,  the  Von  Mises  criterion  is  used  to  define 
failure  in  terms  of  a maximum  allowable  yield  zone  near  the  re-entrant  corner,  with  failure 
occurring  when  the  plate  yields  outside  of  this  zone.  For  plane  strain,  the  radial  stress  at  a 
distance  r from  the  re-entrant  corner  with  angle  0(see  Fig.  4.2),  is  (Williams76  1952) 

a ~ r(-'  (4.1) 

where  is  solved  from  the  following  transcendental  equation: 

Sin(C0)  = ±£Sin{0)  (4.2) 

To  avoid  the  need  for  repetitively  solving  Eq.  (4.2),  a third  order  polynomial  approximation 
of  4"  as  a function  of  9 (in  radians)  was  constructed.  Nineteen  £ values  corresponding  to 
equally  spaced  values  of  9 between  ;rand  2>n/2  were  used  to  obtain: 

£ = 7.457  -4.0460  + 0.806 192  - 0.0549 1<93  (4.3) 

The  %MaxErr  value  (see  Eq.  2.21)  for  the  data  points  used  to  construct  the  approximation  is 
0.41%. 

The  plane  strain  governing  equations,  obtained  from  the  equilibrium  equations,  the 
stress-strain  and  strain-displacement  relationships  (Cook,  Malkus  and  Plesha77  1989,  pp.  20-21), 
depend  on  nine  parameters:  A,  /?,  y,  a,  b,  to,  P,  E and  v.  By  defining  the  non-dimensional 


parameters 
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x = x/a 
u={Eto/P)u 

£*=(Eato/P)£X 

£:=(Eato/P)£: 

r»  = iEat,/p)r„ 


z=z/t0 
v = (£a/P)v 

3s,  = (Mt,/p)ffx 

={Ahto/P)°: 


the  governing  equations  may  be  written  in  non-dimensional  form  as: 

2(1  - v){t0  /a)u,~  +{to  /a)v,~~  +(1  - 2v)(a/t0  )u,n  = 0 

(°A.  )«  ,h  +(1  - 2w)(f.  /fl)v ,_  +2(1  - u)(flr/ro  )v,n  = 0 

The  associated  non-dimensional  boundary  conditions  are  summarized  in  Table  (4.1). 


Table  4.1:  Non-dimensional  boundary  conditions. 


Boundary 

Prescribed  Stress 

x = -y2 

3 = -X  r = 0 

X xz 

x=y2 

za 

X 

n 

« 

II 

O 

z = 0 

Cp 

N 

II 

o 

» 

II 

o 

z - 1 

3,  =0  r„  = 0 

N* 

II 

cp 

N 

II 

O 

II 

o 

(y2  + p)<x  <{y2  + p +y) 

and 

Normal  Stress: 

z = mix  + c 

&N  =° 

m = -(1  - X)l  y 

Surface  Shear: 

\-X  + ip-2Xp  + 2y 
2 y 

o 

II 

* 

(4.4) 


(4.5) 


From  Eq.  (4.5)  and  Table  (4.1),  X,  ft,  y,  a/t0  and  u are  identified  as  the  intrinsic  variables 
associated  with  the  plane  strain  solution.  However,  the  stress  distribution  near  the  re-entrant 
corner  also  depends  on  the  radial  distance  measured  from  the  corner  (see  Eq.  4.1).  An 
investigation  of  several  plate  configurations  revealed  that  the  maximum  stress  always  occurs 
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on  the  top  surface  of  the  thinner  section  of  the  plate  where  crx  is  the  only  non-zero  stress 
component,  and  the  Von  Mises  failure  criterion  reduces  to 

Abt„CTy 


=• 


(4.6) 


where  ay  denotes  the  material  yield  stress.  By  defining  the  non-dimensional  distance 


~ _ r 
t. 


(4.7) 


the  response  surface  approximation  of  the  non-dimensional  stress  distribution  on  the  top 
surface  of  the  thinner  section  of  the  plate  may  be  approximated  as 

ax  = C,  (A,  p,  y)  + C2  (A,  p,  y)r(A  (4.8) 

where  4”  also  depends  on  A,  y and  a/t„  through  the  angle  6.  In  the  present  dissertation, 
response  surface  approximations  were  constructed  to  represent  the  unknown  functions  Ci  and 
C?  of  Eq.  (4.8). 

Thus,  there  are  seven  intrinsic  variables  (A,  /?,  y,  a/t0,  v,  r and  4")  associated  with  the 
solution  for  the  stress  distribution  near  the  re-entrant  corner.  However,  it  is  expected  that  the 
non-dimensional  stress  distribution  varies  only  for  small  values  of  a/t0  (thick  plates),  and  that 
for  large  a/t0  values  (thin  plates),  the  influence  of  a/tQ  may  be  ignored.  Additionally,  it  is  also 
anticipated  that  small  changes  in  Poisson’s  ratio  about  the  baseline  value  of  0.3  should  not 
have  a large  influence  on  the  non-dimensional  stress  distribution.  For  the  case  where  A=0.5, 
P=  0,  y=  0 (corresponding  to  9=7>n/2  and  4=0-55)  and  r =0.1  (10%  of  ta  ),  varying  a/t0 
between  1 and  400  changes  the  non-dimensional  stress  by  less  than  1%  for  a/t0  larger  than  25. 
Furthermore,  for  the  case  where  A=0.5,  /?=0,  y=  0,  r =0.03  and  4=0.55,  varying  the  value  of 
v from  0.2  to  0.4  changes  the  non-dimensional  stress  by  only  1.85%.  In  addition  to  these 
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evaluations  for  a single  design  point,  random  combinations  of  A,  /?  and  y were  tested  as 
discussed  subsequently  in  Section  4.3. 1.1. 

Hence,  when  t>«0.3  and  a/to  >25,  the  Sx  stress  on  the  top  surface  of  the  thin 
section  of  the  plate  may  be  described  in  terms  of  only  four  predictor  variables  A,  fi,  yand  rc~  , 
where  r and  4” are  combined  into  a single  variable  riA  . Recall  that  predictor  variables  are  the 
variables  of  the  response  surface  approximation  used  to  approximate  &x  . 


4.2.2. Buckling  Load  of  the  Plate 

The  buckling  load  of  the  plate  depends  on  eight  parameters:  A,  J3,  y,  a,  b,  t0,  E and  u. 
The  flexural  stiffness  of  the  plate  is  given  by 


D = 


Et 3 

12(1 -v2) 


(4.9) 


where  the  thickness  t is  a function  of  t0,  A,  ft,  r and  x.  The  non-dimensional  parameters 


are  defined,  and 


~ x 
x — — 

Q 


y = 


y_ 

a 


7-  1 

a 


w 


w = - 


N = 


N. 


Euler 


d=d 


Ea 3 12(1 -t>2) 


(4.10) 


^ Euler 


7t2D 


a 


(4.11) 


where  D is  a representative  value  of  D,  defined  subsequently  in  Eq.  (4.19).  Then  the 
governing  equation  and  boundary  conditions  for  the  buckling  load  due  to  uniform  axial 
compression  (Timoshenko78  1936,  p.  297  and  pp.  337-350)  may  be  non-dimensionalized  to 
obtain 
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(Dw~  L + 2(^W  l + fa,*  U * L “ "fear  L 

+ 2u(5v?  ~ )_  = ~n2DN xw~  (4.12) 

where: 

For  x = -o/2  and  x=a/2:  w = 0 w~=0 

b b 

F°r  ^ = ~~  2^  and  ^ = 2^ : + ^ = ° ^ + <2  ~ ”)*«?  = 0 (4-13) 

The  value  of  Nx  corresponding  to  the  buckling  load  is  denoted  by  NCril.  Five 
intrinsic  variables  associated  with  the  solution  for  the  buckling  load  are  identified  as:  v,  a/b,  A, 
(J  and  y.  The  non-dimensionalization  of  Eq.  (4.10)  is  based  on  Euler  buckling  of  the  plate. 
The  non-dimensional  buckling  load  is  independent  of  Poisson’s  ratio  and  a/b  because  the 
Euler  buckling  load  is  independent  of  these  parameters.  However,  for  large  values  of  a/b  the 
structure  behaves  like  a beam,  rather  than  a plate,  with  the  buckling  load  being  higher  by  a 
factor  of  l/(l-t>2) . The  dependence  of  the  non-dimensional  buckling  load  on  the  variables  o 
and  a/b  was  tested  by  varying  these  variables  for  the  case  2,=0.5,  J3=0  and  y =0.  It  was  found 
that  NCrjl  changes  by  3.90%  when  v is  varied  from  0.2  to  0.4  (a/b= 1),  while  NCrit  changes  by 
7.5%  when  a/b  is  varied  from  0.5  and  10  (v—0.3).  Note  that  the  relatively  large  error  of  7.5% 
is  due  to  the  fact  that  the  structure  behaves  like  a beam  when  a/b= 10.  Additional  evaluations 
of  the  magnitude  of  the  dependence  of  the  buckling  load  on  Poisson’s  ratio  and  a/b  at  other 
random  design  points  are  described  subsequently  in  Section  4.3. 1.2. 

Hence,  for  a plate  with  a Poisson’s  ratio  value  close  to  0.3,  the  response  surface 
approximation  for  NCrjl  has  only  three  predictor  variables  2, 1 3 and  y. 
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4.3.  Response  Surface  Approximation  Construction 
The  objective  of  the  current  design  problem  is  to  maximize  a safety  measure  of  the 
plate  with  change  in  thickness  for  a given  weight.  In  the  present  dissertation,  the  factor  of 
safety  is  the  safety  measure  considered  for  deterministic  design,  while  one  minus  the 
possibility  of  failure  is  the  safety  measure  considered  for  fuzzy  set  based  designing  for 
uncertainty.  When  considering  both  the  yielding  (Von  Mises  yield  stress  failure  criterion) 
failure  mechanism  and  the  buckling  response  constraint,  the  failure  load  Pf  of  the  plate  is 
given  by 


Pf  = Min{ 


Yield  Stress 
[Buckling  Load 


with  failure  defined  to  occur  when: 


P-Pf>  0 


(4.14) 


(4.15) 


In  the  present  example  problem,  response  surface  approximations  form  an  integral 
part  of  the  fuzzy  set  based  design.  As  outlined  in  Fig.  (4.5),  two  levels  of  response  surface 
approximations  are  employed  during  different  stages  of  the  design  process.  On  the  first  level, 


First  Level  Response 
Surface  Approximation: 

Failure  Criteria 


Second  Level  Response 
Surface  Approximation: 

Possibility  of  Failure 


Figure  4.5: 


Summary  of  two  level  response  surface  approximation  construction. 
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response  surface  approximations  of  the  non-dimensional  stress  distribution  crx  and  non- 
dimensional  buckling  load  NCrjl  are  constructed  and  are  used  to  replace  computationally 
expensive  finite  element  analyses  during  the  evaluation  of  the  possibility  of  failure  of  the 
structure.  These  response  surface  approximations  are  constructed  based  on  finite  element 
analyses  and  are  used  mainly  to  reduce  the  computational  burden  associated  with  calculating 
the  possibility  of  failure  using  the  vertex  method. 

On  the  second  level,  response  surface  approximations  of  the  possibility  of  failure  as  a 
function  of  the  nominal  values  of  the  design  variables  and  the  level  of  uncertainty  associated 
with  these  variables  are  constructed.  These  approximations  are  constructed  based  on  results 
obtained  from  the  vertex  method,  using  the  first  level  response  surface  approximations.  The 
second  level  response  surface  approximations  are  constructed  mainly  to  simplify  the 
integration  of  the  analysis  code  with  the  optimization  algorithm,  as  well  as  to  filter  out 
numerical  noise  by  providing  a smooth  approximate  response  function.  Filtering  out  noise  in 
the  approximate  response  function  allows  the  use  of  a derivative-based  optimization 
algorithm.  The  generalized  reduced  gradient  algorithm  provided  with  Microsoft  Excel 
Version  8.069  (1997)  was  used  in  the  present  example  problem. 

4.3.1. Stress  Distribution  and  Buckling  Load  Response  Surface  Approximations 

From  Section  4.2.1,  the  stress  distribution  response  surface  approximation  may  be 

written  in  functional  form  as  follows: 

cr,  =crx(A,/},y,r(-')  (4.16) 

The  behavior  of  the  stress  distribution  near  the  re-entrant  corner  of  the  plate  was 
determined  by  varying  each  of  the  four  predictor  variables  individually  for  specific 
combinations  of  the  other  three  variables.  It  was  determined  that  a fourth  order  polynomial 
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in  the  four  predictor  variables  (70  parameters)  was  suitable  for  an  initial  response  surface 
approximation. 

The  buckling  load  response  was  similarly  investigated  and  the  resulting  graphs 
indicated  that  the  behavior  of  the  buckling  load  response  is  more  complex  than  the  stress 
distribution  response.  The  additional  complexity  is  mainly  due  to  the  presence  of  two 
buckling  modes.  A local  buckling  mode  is  active  for  high  slenderness  ratios  of  the  thin  section 
of  the  plate,  while  a global  buckling  mode  is  active  otherwise.  To  model  these  two  buckling 
modes  accurately,  it  was  necessary  to  construct  two  independent  response  surface 
approximations,  one  for  each  buckling  mode.  These  response  surface  approximations  may  be 
written  in  functional  form  (from  Section  4.2.2)  as: 

N Local  = NLocal{h,P,Y)  ^Global  ~ ^GIobal(^’ P ’T') 

The  necessity  of  two  separate  response  surface  approximations  implies  that  the  buckling  load 
data,  obtained  from  the  finite  element  analyses,  must  be  divided  into  two  groups 
corresponding  to  the  two  buckling  modes.  A simple  geometric  criterion,  which  assigns  data 
points  to  the  appropriate  buckling  mode  response  surface  approximation,  was  defined  based 
on  the  slenderness  ratio  of  the  thin  section  of  the  plate.  This  criterion  is  written  in  non- 
dimensional  form  as: 


Buckling  mode  = 


Local  mode  if  (0.5  /?  y)  > q ^ 
X 

(0.5 -P-y) 

Global  mode  if  , <0.6 


(4.18) 


Two  separate  definitions  of  the  D value  of  the  Euler  buckling  load  used  in  the  non- 
dimensionalization  of  Eq.  (4.10)  was  defined  based  on  the  local  and  global  buckling  modes 
respectively.  For  the  local  buckling  mode,  the  thickness  of  the  thin  section  of  the  plate,  Xta , 
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was  used  as  an  equivalent  thickness  value,  while  the  average  t 3 value  was  used  as  an  equivalent 
thickness  value  for  the  global  buckling  mode.  The  D value  is  then  defined  as: 


D = 


— if  local  buckling  mode 
12(1  -v2) 

12(1  ^ §1°^  buckling  mode 


(4.19) 


After  investigating  several  different  order  polynomials,  it  was  determined  that  a third 
order  polynomial  (20  parameters)  is  suitable  as  a response  surface  approximation  for  the  local 
buckling  mode,  and  a fourth  order  polynomial  (35  parameters)  is  suitable  as  a response  surface 
approximation  for  the  global  buckling  mode. 

The  design  space  used  for  constructing  the  stress  distribution  and  buckling  load 
response  surface  approximations  is  summarized  in  Table  (4.2)  and  is  shown  graphically  for  X , 
/land  yin  Fig  (4.6).  In  Table  (4.2)  the  upper  limit  on  r limits  the  radius  of  the  yield  zone 


Table  4.2:  Design  space  for  constructing  the  response  surface  approximations. 


Response  Variable  Range 

^ 0.2  < T<  1.0 


- 0.475  <p<  0.475 
0 < y < 0.475  -/? 


Figure  4.6: 


Three  dimensional  representation  of  the  design  space. 
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about  the  re-entrant  corner  to  be  no  greater  than  80%  of  the  thickness  of  the  thin  section  of 
the  plate.  The  upper  bound  on  y is  dictated  by  the  geometry  of  the  transition  region. 

For  the  isotropic  plate  example  problem,  each  finite  element  analysis  was 
computationally  inexpensive  to  evaluate.  Since  it  was  feasible  to  perform  a large  number  of 
numerical  experiments,  it  was  not  necessary  to  use  statistical  design  of  experiments  for 
selecting  a small  number  of  data  points  at  which  to  evaluate  the  response  function.  Instead, 
the  response  function  was  evaluated  at  a large  number  of  data  points.  Although  this  approach 
is  rarely  practical,  it  provided  a large  number  of  data  points  that  were  used  as  an  independent 
data  set  for  evaluating  the  predictive  capabilities  of  the  response  surface  approximations.  This 
independent  data  set  was  used  to  validate  the  proposed  methodology  before  considering  the 
second,  more  complicated,  example  problem.  The  numerical  experiments  were  conducted  for 
a three-dimensional  grid  of  equally  spaced  data  points  with  eight  data  points  in  each  of  the  A,  /? 
and  y directions.  For  a rectangular  domain,  this  procedure  would  yield  a total  of  512  data 
points.  However,  due  to  the  prism  or  wedge-like  domain  of  Fig.  (4.6),  only  288  of  these  512 
data  points  fell  within  the  current  design  space. 

4.3.1. 1.  Stress  distribution  response  surface  approximation 

For  the  stress  distribution,  each  of  the  288  plate  configurations  (corresponding  to  288 

finite  element  analyses)  included  a number  of  data  points  with  different  rh  1 values 
(corresponding  to  different  finite  elements).  The  stress  values  for  all  finite  elements  on  the  top 
surface  of  the  thin  section  of  the  plate  were  recorded  for  node  points  within  80%  of  the 
thickness  of  the  thin  section  of  the  plate,  measured  from  the  re-entrant  corner.  This  process 
yielded  a total  of  3,023  data  points  that  were  divided  into  two  groups  by  assigning  a random 
number  between  0 and  1 to  each  data  point.  Data  points  with  a random  number  less  than  0.7 


56 


were  allocated  to  the  set  of  data  points  used  to  construct  the  response  surface  approximation, 
yielding  2,124  data  points  (Data  Set  1).  The  remaining  899  data  points  were  reserved  to 
evaluate  the  predictive  capabilities  of  the  response  surface  approximation  (Data  Set  2). 

Using  the  mixed  stepwise  regression  procedure,  the  optimum  number  of  parameters  in 
a reduced  response  surface  approximation  for  the  stress  distribution  was  found  to  be  43. 
Using  both  Data  Set  1 and  Data  Set  2,  the  predictive  capabilities  of  the  initial  response  surface 
approximation  (70  parameters)  and  of  the  reduced  response  surface  approximation 
(43  parameters)  can  be  compared  by  comparing  the  results  shown  in  Table  (4.3). 


Table  4.3:  Stress  distribution  response  surface  approximation. 


Model 

Cp 

R2 

Adj-R2 

%RMSE 

%RMSE 

PRESS 

%Avg 

Err 

%RMSE 

%Max 

Err 

Stress 

Data  Set  1 (2,124  Data  Points) 

Data  Set  2 (899  Data  Points) 

Full 

(70  terms) 

70 

0.9984 

0.9983 

3.2478 

3.4172 

2.4552 

3.5934 

16.443 

Reduced 
(43  terms) 

31.927 

0.9983 

0.9982 

3.2964 

3.3886 

2.4339 

3.4754 

17.120 

The  results  in  Table  (4.3)  indicate  that  both  models  have  similar  predictive  capabilities. 
The  reduced  response  surface  approximation  is  slightly  more  accurate  (except  for  the 
%MaxErr  values),  and  is  preferred  since  it  represents  a smoother  function  due  to  a smaller 
number  of  parameters.  The  error  terms  obtained  from  Data  Set  1 are  slightly  optimistic  due 
to  the  fact  that  they  are  calculated  from  data  points  used  to  construct  the  response  surface 
approximation.  However,  the  results  obtained  from  the  %RMSEpress  value  compare  well  with 
the  results  obtained  from  Data  Set  2.  This  trend  is  a general  trend  when  the  number  of  data 
points  is  much  larger  than  the  number  of  parameters  in  the  response  surface  approximation. 
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Figure  (4.7)  shows  the  correlation  between  the  non-dimensional  finite  element  stress 
values  and  the  non-dimensional  stress  values  obtained  from  the  response  surface 
approximation.  The  correlation  for  the  data  points  of  Data  Sets  1 and  2 is  shown  as  circular 
symbols.  If  there  were  100%  correlation  between  the  finite  element  and  the  response  surface 
approximation  stress  values,  the  circular  symbols  would  lay  on  the  ideal  diagonal  line.  From 
Fig.  (4.7)  it  is  clear  that  the  circular  symbols  are  close  to  the  ideal  diagonal  line,  indicating  a 
high  correlation  or  accuracy  between  the  finite  element  and  response  surface  approximation 
stress  values. 


Figure  4.7:  Response  surface  approximation  compared  to  finite  element  non-dimensional 

stress  values  (Data  Sets  1 and  2). 

The  ability  of  the  response  surface  approximation  to  predict  accurate  values  for 
Poisson’s  ratios  that  are  different  from  0.3  was  evaluated  with  ten  random  combinations  of  A, 
/3  and  y.  The  values  u=0.25  and  o=0.35  were  used  for  each  configuration.  The  average  error 
for  this  evaluation  set  was  3.3%,  indicating  that  small  changes  in  Poisson’s  ratio  do  no  effect 
the  results. 

An  important  and  desirable  characteristic  of  response  surface  approximations  is  their 
tendency  to  filter  out  numerical  noise.  This  characteristic  is  illustrated  for  the  maximum 
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stress  value  on  the  top  surface  of  the  thin  section  of  the  plate  (independent  of  r ) as  a function 
of  X for  P=y=  0.  These  values  of  /3  and  /correspond  to  a sharp  transition  at  the  mid-length  of 
the  plate.  For  this  illustration,  60  data  points  were  evaluated  and  the  experimental  results  are 
shown  as  circular  symbols  in  Fig.  (4.8).  The  response  obtained  from  the  response  surface 
approximation  is  shown  as  a solid  line  in  Fig.  (4.8). 


Figure  4.8:  Maximum  non-dimensional  stress  value  as  a function  of  X (/?—  0,  y=  0). 

There  are  some  inconsistencies  in  the  data  obtained  from  the  finite  element  model, 
which  may  be  attributed  to  the  mesh  generator  that  did  not  provide  finite  element  models 
with  a consistent  level  of  fidelity.  In  particular,  the  mesh  generator  yielded  less  accurate  finite 
element  models  when  y is  equal  to  0 and  X is  close  to  1,  corresponding  to  a sharp  transition  in 
thickness  with  a small  change  in  thickness.  The  results  presented  in  Fig.  (4.8)  indicate  that  the 
response  surface  approximation  eliminated  these  inconsistencies,  or  noise,  and  indicate  that 
not  all  differences  between  predicted  and  experimental  response  values  should  be  considered 
to  be  errors.  For  some  cases,  the  predicted  response  values  may  indeed  be  more  accurate  than 
the  experimental  response  values  as  is  the  case  in  Fig.  (4.8)  when  X is  close  to  1.  When  X is 
close  to  1,  the  non-dimensional  finite  element  stress  values  should  approach  a value  equal  to  1 
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and  are  thus  in  error.  The  problem  with  the  mesh  generator  disappears  when  X is  equal  to  1 
(no  change  in  thickness)  since  the  non-dimensional  finite  element  stress  value  for  this  design 
point  is  also  equal  to  1,  corresponding  to  no  stress  concentration. 

4.3. 1.2.  Buckling  load  response  surface  approximation 

For  representing  the  buckling  load,  Eq.  (4.18)  was  used  to  divide  the  288  data  points 
into  126  data  points  for  the  construction  of  the  local  buckling  mode  response  surface 
approximation,  and  162  data  points  for  the  global  buckling  mode  response  surface 
approximation.  An  additional  set  of  196  data  points  was  generated  to  be  equally  spaced, 
midway  between  the  existing  288  data  points,  with  83  of  these  data  points  used  for  evaluating 
the  local  buckling  mode  response  surface  approximation,  and  113  of  these  data  points  used  for 
evaluating  the  global  buckling  mode  response  surface  approximation. 

Using  the  mixed  stepwise  regression  procedure,  the  best  reduced  response  surface 
approximation  for  the  local  buckling  load  (cubic  polynomial)  has  19  parameters  while  the 
reduced  response  surface  approximation  for  the  global  buckling  load  (quartic  polynomial)  has 
25  parameters.  The  predictive  capabilities  of  the  models  are  summarized  in  Table  (4.4),  and  a 
graphical  representation  of  the  results  is  shown  in  Fig.  (4.9).  Figure  (4.9)  shows  the 
correlation  between  the  finite  element  and  response  surface  approximation  non-dimensional 
buckling  load  values  for  the  data  points  of  Data  Sets  1 and  2.  As  was  the  case  in  Fig.  (4.7),  the 
high  accuracy  of  the  non-dimensional  buckling  load  response  surface  approximations  is 
illustrated  by  the  fact  that  the  symbols  representing  the  actual  data  points  are  very  close  to  the 
ideal  diagonal  line. 

In  addition,  to  determine  the  accuracy  of  the  response  surface  approximations  for 
different  values  of  o and  a/b,  ten  random  combinations  of  X,  f3  and  y were  tested  for  t>=0.25 
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Table  4.4:  Buckling  modes  response  surface  approximations. 


Model 

Cp 

R2 

Adj-R2 

%RMSE 

%RMSE 

PRESS 

%RMSE 

Err 

%Max 

Err 

Local 

Buckling 

Data  Set  1 (126  Data  Points) 

Data  Set  2 (83  Data  Points) 

Full 

(20  terms) 

20 

0.9999 

0.9998 

0.5541 

0.6927 

0.7056  0.9567 

2.8494 

Reduced 
(19  terms) 

19.402 

0.9999 

0.9998 

0.5550 

0.6920 

0.7055  0.9598 

2.8792 

Global 

Buckling 

Data  Set  1 (162  Data  Points) 

Data  Set  2 (113  Data  Points) 

Full 

(35  terms) 

35 

0.9912 

0.9889 

2.5538 

3.1695 

1.1998  1.6895 

5.4340 

Reduced 
(25  terms) 

18.122 

0.9910 

0.9895 

2.4888 

3.0202 

1.1936  1.6935 

5.8767 

Figure  4.9:  Response  surface  approximation  compared  to  finite  element  non-dimensional 

buckling  loads  for  both  buckling  modes  (Data  Sets  1 and  2). 

and  t>=0.35,  and  for  a/b= 0.1  and  a/b- 5.  For  a/b=0.\  the  average  error  was  4.5%  while  for 
a/b= 5 the  average  error  ranged  from  2.7%  for  u—0.25  to  9.1%  for  o=0.35.  This  larger  error 
reflects  the  fact  that  the  structure  behaves  like  a beam  rather  than  a plate  for  a/b= 5. 

As  was  the  case  for  the  stress  distribution  response  surface  approximation,  both  the 
reduced  and  the  full  response  surface  approximations  provide  similar  predictive  capabilities. 
As  before,  the  reduced  response  surface  approximations  are  preferred  over  the  full  response 
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surface  approximations,  since  the  reduced  response  surface  approximations  represent 
smoother  functions  due  to  a smaller  number  of  parameters. 

4.3.2. Possibility  of  Failure  Response  Surface  Approximations 

The  A,  P and  /design  space  of  Table  (4.2)  was  also  used  to  construct  the  possibility  of 
failure  response  surface  approximation,  with  numerical  experiments  conducted  on  an  evenly 
spaced  grid  consisting  of  11  data  points  in  each  of  the  A,  P and  / directions.  Additionally, 
seven  evenly  spaced  levels  of  uncertainty  were  considered  between  ± 2%  and  ± 20%,  yielding  a 
total  of  2,629  data  points  in  the  design  space.  At  each  data  point  the  possibility  of  failure 
values  due  to  yielding  and  due  to  violating  the  buckling  load  constraint  were  evaluated  from 
Eqs.  (4.14)  and  (4.15)  using  the  vertex  method  and  the  previously  developed  reduced  stress 
distribution  and  buckling  load  response  surface  approximations. 

Two  response  surface  approximations,  according  to  the  two  possibilities  of  failure, 
were  constructed  using  all  data  points  with  possibility  of  failure  not  equal  to  either  0 or  1. 
This  resulted  in  499  data  points  for  constructing  the  yield  stress  failure  mode  response  surface 
approximation  and  573  data  points  for  the  buckling  load  constraint  failure  response  surface 
approximation.  The  resulting  predicted  possibility  of  failure  is  then  obtained  from 

n(p-p f){A,p,y,u)  = nun  (id  TieUstna  (A,  p,  y,u),fl  Bucklmg  (2,  /?,  /,  a ))  (4.20) 

where  u denotes  the  uncertainty  associated  with  the  design  variables  and,  as  before,  the  caret 
symbol  refers  to  predicted  values. 

It  was  found  that  a general  fourth  order  polynomial  (70  parameters)  resulted  in 
accurate  response  surface  approximations  for  both  possibility  of  failure  response  surface 
approximations.  These  general  response  surface  approximations  were  reduced  using  the 
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mixed  stepwise  regression  procedure  and  the  Cp  statistic.  The  predictive  capabilities  of  the 
response  surface  approximations  are  summarized  in  Table  (4.5). 

Table  4.5:  Possibility  of  failure  response  surface  approximation. 

Model  Cp  R2  Adj-R2  RMSE  PRESS 
Stress  Fourth  Order  Model  (499  Data  Points) 

70.000  0.9988  0.9986  2.2525  2.6205 

53.027  0.9988  0.9986  2.2371  2.5205 

Fourth  Order  Model  (573  Data  Points) 

70.000  0.9982  0.9980  2.7342  3.0991 

48.649  0.9982  0.9980  2.7118  3.0211 


For  the  possibility  of  failure  response  surface  approximations,  only  a single  data  set 
was  used  to  construct  the  approximations.  In  this  case,  the  %RMSEPress  value  was  used  to 
compensate  for  the  possibly  overly  optimistic  error  estimates.  As  noted  previously 
(see  Section  4.3. 1.1),  the  %RMSEPREss  value  provides  good  results  when  the  number  of  data 
points  used  to  construct  the  response  surface  approximation  is  much  larger  than  the  number 
of  parameters  in  the  response  surface  approximation.  Since  only  one  data  set  was  used  to 
construct  the  response  surface  approximations,  the  predictive  capabilities  of  the  response 
surface  approximations  are  not  presented  in  graphical  form. 


Full 

(70  Terms) 
Reduced 
(57  Terms) 


Full 

(70  Terms) 
Reduced 
(59  Terms) 

Buckling 


CHAPTER  5 

ISOTROPIC  PLATE  OPTIMIZATION 

The  design  for  uncertainty  problem  considered  here,  involves  maximizing  a safety 
measure  (i.e.,  the  factor  of  safety  for  deterministic  design  and  one  minus  the  possibility  of 
failure  for  fuzzy  set  based  design)  of  an  aircraft  structure  for  a given  weight.  The  structure 
considered  is  typical  of  a wing  skin  application  where  the  thickness  is  reduced  from  the  root 
to  the  tip,  and  is  subject  to  uncertainty  in  the  geometric  variables,  the  material  properties  and 
the  applied  load.  It  is  assumed  that  the  structure  will  be  built  well  into  the  future,  from 
materials  not  yet  available.  In  this  case  little  information  on  the  uncertainty  associated  with 
the  uncertain  problem  parameters  is  available,  and  the  uncertainty  is  based  on  expert  opinion 
and  assumptions  made  by  the  designer.  Fuzzy  set  theory  is  appropriate  to  model  uncertainty 
with  little  available  information  and  is  used  in  the  present  dissertation.  Response  surface 
approximations  are  used  throughout  the  design  process,  mainly  to  reduce  the  computational 
burden  associated  with  designing  for  uncertainty,  but  also  to  integrate  the  analysis  code  with 
the  optimization  algorithm.  Additionally,  response  surface  approximations  filter  out 
numerical  noise  in  the  approximate  response  function  used  during  the  optimization  process, 
thus  allowing  the  use  of  a derivative-based  optimization  algorithm. 

5.1.  General  Design  Problem 

An  isotropic  plate  with  a change  in  thickness  across  its  width  (see  Fig.  5.1)  is  the 
present  design  problem.  All  problem  parameters  are  assumed  to  be  uncertain.  Yielding  is 
considered  as  a failure  mechanism,  while  the  plate  is  also  designed  to  be  buckling  resistant. 
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Figure  5.1:  Cross-section  model  with  non-dimensional  variables  shown. 


The  Von  Mises  failure  criterion  is  used  to  calculate  the  failure  load  of  the  plate  due  to 
yielding.  For  the  buckling  load  constraint,  the  design  fails  if  the  applied  load  exceeds  the 
buckling  load  constraint.  Based  on  these  two  criteria,  the  failure  load  Pf  is  calculated  from 


Pf  = Min< 


°rMto 
_2  77 l/  \3 


NCrl,n2Eb(AtoY 


12(1  -u2)a2 


with  failure  defined  to  occur  when: 


(5.1) 


P-Pf>  0 (5.2) 

Note  that  throughout  the  design  problem  the  buckling  load  Nat  is  non-dimensionalized  with 
respect  to  the  Euler  buckling  load  of  a plate  with  a constant  thickness  equal  to  X to . 

The  objective  of  the  design  problem  is  to  maximize  a safety  measure  of  the  structure 
for  a given  weight.  The  results  obtained  from  a traditional  deterministic  approach,  using  a 
factor  of  safety  to  account  for  the  uncertainty,  are  compared  to  the  results  obtained  from  a 
fuzzy  set  based  approach.  For  the  deterministic  design,  the  factor  of  safety  was  used  as  the 
safety  measure,  while  one  minus  the  possibility  of  failure  was  used  as  the  safety  measure  for 
the  fuzzy  set  based  design.  In  addition,  the  fuzzy  set  based  analysis  was  also  used  to  study  the 
dependence  of  the  weight  of  the  final  design  on  the  level  of  uncertainty  associated  with  the 
design  variables  X,  /?  and  y,  with  the  results  presented  in  the  form  of  a cost  versus  weight 
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design  chart.  Different  levels  of  uncertainty  for  the  design  variables  were  considered,  since 
these  geometric  variables  have  the  largest  influence  on  the  manufacturing  cost  of  the  plate.  If 
the  tolerances  of  these  variables  can  be  relaxed  without  a large  weight  penalty,  substantial  cost 
reductions  can  be  achieved  in  manufacturing  the  plate.  The  problem  parameters  and 
associated  levels  of  uncertainty  used  in  the  example  problem  are  summarized  in  Table  (5.1). 


Table  5.1:  Uncertain  problem  parameters. 


Variable 

Nominal  Values 

Level  of 
Uncertainty,  u 

k* 

[0.2 -1.0] 

±5% 

P 

o 

1 

-r 

o 

±5% 

r ' 

[0  - 0.8] 

±5% 

a 

90  in 

±5% 

b 

50  in 

±5% 

to 

3 in 

±5% 

E 

30  000  ksi 

±5% 

V 

0.29 

±5% 

Oy 

26  ksi 

±10% 

r 

5 A.  ta 

±10% 

P 

72  5 000  lb 

±10% 

t Design  Variables 


It  is  important  to  note  that  there  exists  fundamental  differences  in  the  manner  the 
deterministic  and  fuzzy  set  based  approaches  maximize  the  safety  measure  for  a given  weight. 
The  deterministic  design  tends  to  equalize  the  two  failure  loads  of  Eq.  (5.1),  while  the  fuzzy 
set  based  design  tends  to  equalize  the  possibility  of  failure  values  due  to  yielding  and  due  to 
violating  the  buckling  load  constraint.  If  the  design  variable  limits  are  not  active  at  the 
optimum  design,  one  would  expect  a deterministic  design  with  equal  failure  load  values  and  a 
fuzzy  set  based  design  with  equal  possibility  of  failure  values  for  the  two  criteria  considered  in 
the  design. 
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5.2.  Deterministic  Design  Problem 

In  the  deterministic  design,  a factor  of  safety  is  used  to  account  for  the  uncertainty  in 
the  problem  parameters.  The  safety  measure  of  the  structure  is  thus  maximized  by 
maximizing  the  factor  of  safety  for  a given  weight.  For  the  deterministic  design,  failure  is 
calculated  from  Eqs.  (5.1)  and  (5.2),  using  the  nominal  values  of  the  problem  parameters,  while 
the  non-dimensional  stress  distribution  and  buckling  load  values  are  obtained  from  the 
reduced  response  surface  approximations  of  CHAPTER  4. 


5.2.1. Optimization  Problem 

The  objective  of  the  design  problem  is  to  maximize  the  factor  of  safety  for  a given 
weight.  However,  since  it  is  difficult  to  specify  a meaningful  upper  limit  of  the  weight,  it  was 
decided  to  minimize  the  weight  of  the  plate,  for  a given  factor  of  safety,  in  the  deterministic 
design.  The  resulting  minimum  weight  was  then  used  as  the  upper  limit  of  the  weight  for  the 
fuzzy  set  based  design.  It  was  assumed  that  a factor  of  safety  of  1.5  would  be  adequate  to 
account  for  the  uncertainty  associated  with  the  problem  parameters  as  summarized  in 
Table  (5.1).  The  non-dimensional  cross-sectional  area  A of  the  plate  was  used  as  a 
representative  value  of  the  weight  and  the  resulting  optimization  problem  was  defined  as: 
Minimize: 


Subject  to: 


a =U\+2/3+r)+^(i-2/3-r) 

At  „ Z z 


(5.3) 


A 

0.4 


+ 1>0 


A 

0.4 


>0 


y > 0 


0.4 


Pf 

— -1.5  > 0 
p 


(5.4) 
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The  constraints  involving  /?  and  y are  geometric  constraints  and  Pf  is  calculated  from  Eq.  (5.1) 
using  the  nominal  values  of  the  design  variables  and  the  reduced  stress  distribution  and 
buckling  load  response  surface  approximations. 


5.2.2. Results 

The  minimum  non-dimensional  cross-sectional  area  of  the  plate,  obtained  from 
Eqs.  (5.3)  and  (5.4),  are  summarized  in  Table  (5.2).  For  the  deterministic  optimum  design, 
both  the  Von  Mises  yield  stress  failure  criterion  and  the  bucking  constraint  are  active  at  the 
optimum  design  point.  The  optimum  design  corresponds  to  a plate  with  a change  in  thickness 
that  starts  at  the  minimum  allowable  distance  from  the  left  endpoint  of  the  plate  (see  Fig.  5.1) 
with  a very  short  transition  zone  (small  y value).  Both  the  possibility  of  failure  values 
calculated  from  the  vertex  method  and  those  predicted  by  the  reduced  possibility  of  failure 
response  surface  approximations  (values  in  parenthesis)  are  summarized  in  Table  (5.2).  It  is 
clear  that  even  though  both  the  Von  Mises  yield  stress  failure  criterion  and  the  buckling  load 


Table  5.2:  Deterministic  optimum  (±5%  uncertainty  in  the  design  variables). 


Variable 

Value 

X 

0.6287 

P 

-0.4000 

7 

0.0447 

A' 

0.6741 

Factor  of  safety 

1.5 

(Yield  Stress) 

Factor  of  safety 
(Buckling  Load) 

1.5 

^Yields  tress 

0.0977 

(0.1181) 

^Buckling 

0.3411 

(0.3331) 

Values  in  parenthesis  are  predicted  by  possibility  of  failure  response  surface  approximation. 
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constraint  are  active  at  the  optimum  design  point,  there  exists  a large  difference  between  the 
possibility  of  failure  values  for  these  two  criteria,  with  the  buckling  load  constraint  being 
critical.  The  accuracy  of  the  possibility  of  failure  response  surface  approximation  is 
demonstrated,  since  the  difference  between  the  critical  predicted  and  calculated  (using  the 
vertex  method)  possibility  of  failure  values  at  the  optimum  design  is  only  2.3%.  The 
difference  between  the  predicted  and  calculated  possibility  of  failure  values  according  to  the 
Von  Mises  yield  stress  criterion  is  equal  to  20.9%,  but  yielding  is  not  the  active  failure  mode. 

5.3.  Fuzzy  Set  Based  Design  Problem 

The  fuzzy  set  based  design  problem  uses  fuzzy  set  theory  to  model  the  uncertainty  in 
the  problem  parameters.  The  safety  measure  of  the  structure  is  maximized  by  minimizing  the 
possibility  of  failure,  using  the  optimum  non-dimensional  cross-sectional  area  obtained  from 
Eqs.  (5.3)  and  (5.4)  and  summarized  in  Table  (5.2)  as  an  upper  limit  of  the  weight.  In  the 
fuzzy  set  based  design  problem,  the  possibility  of  failure  is  obtained  from  the  reduced 
possibility  of  failure  response  surface  approximations  of  CHAPTER  4. 

5.3.1. Optimization  Problem 

The  optimization  problem  for  minimizing  the  possibility  of  failure  subject  to  a given 
weight,  may  be  written  as: 

Minimize: 

H^-py)  =n(M/)(A,Ay)  (5.5) 


^-+1>0 

0.4 


Subject  to: 


0.4 
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y > 0 


0.4 


A(A,j3,y) 


-1  = 0 


(5.6) 


Where  niP_p  , denotes  the  possibility  of  failure,  ri(p_p  , denotes  the  approximate 


possibility  of  failure  obtained  from  the  reduced  response  surface  approximations  and  A" 
denotes  the  optimum  non-dimensional  cross-sectional  area  obtained  from  the  deterministic 
design.  Additionally,  the  variables  A,  fi  and  y denote  the  nominal  values  of  the  design 
variables. 


5.3.2. Results 

The  equivalent  fuzzy  set  based  design  using  the  A'  value  of  Table  (5.2)  is  summarized 
in  Table  (5.3).  The  luzzy  set  based  optimum  design  corresponds  to  a plate  where  the  change 
in  thickness  starts  at  the  minimum  allowable  distance  from  the  left  endpoint  of  the  plate  with 
no  transition  zone  (y  value  equal  to  0).  The  fuzzy  set  based  design  eliminates  the  weight  of 


Table  5.3:  Fuzzy  set  based  optimum  (±5%  uncertainty  in  the  design  variables). 


Variable 

Value 

A 

0.6379 

P 

-0.4000 

r 

0.0000 

A’ 

0.6741 

Factor  of  safety 

1.4898 

(Yield  Stress) 

Factor  of  safety 
(Buckling  Load) 

1.5565 

^Yields  tress 

0.1319 

(0.1262) 

^Buckling 

0.2788 

(0.2721) 

Values  in  parenthesis  are  predicted  by  possibility  of  failure  response  surface  approximation. 
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the  ramp  and  uses  it  to  thicken  the  thin  section  of  the  plate.  The  result  is  an  increase  in  the 
stress  concentration  and  an  improvement  in  the  buckling  load  of  the  plate.  The  fuzzy  set 
based  design  thus  attempts  to  equalize  the  possibility  of  failure  values  due  to  yielding  and  due 
to  violating  the  buckling  load  constraint  by  increasing  both  the  stress  concentration  and  the 
buckling  load  of  the  plate.  However,  in  the  present  design  the  design  variable  limits  kept  the 
possibility  of  failure  values  from  becoming  equal  in  the  final  design.  As  before,  both  the 
possibility  of  failure  values  calculated  from  the  vertex  method  and  those  predicted  from  the 
reduced  models  of  CHAPTER  4 (values  in  parentheses)  are  shown  in  Table  (5.3).  The 
difference  between  the  critical  predicted  and  calculated  possibility  of  failure  values  of  the  plate 
at  the  optimum  design  is  equal  to  2.4%.  The  difference  between  the  predicted  and  calculated 
possibility  of  failure  values  according  to  the  Von  Mises  yield  stress  failure  criterion  is  equal  to 
4.3%,  but  again  yielding  is  not  the  active  failure  mode. 

The  membership  functions  of  the  possibilities  of  failure  according  to  the  fuzzy 
function  (P  - P/ ) are  shown  graphically  in  Fig.  (5.2)  for  the  optimum  designs  obtained  from 
the  deterministic  and  fuzzy  set  based  approaches.  The  fuzzy  function  (P  - P/)  is  calculated 
from  Eqs.  (5.1)  and  (5.2),  corresponding  to  the  Von  Mises  yield  stress  and  the  buckling  load 
constraint  criteria.  The  membership  functions  of  the  optimum  deterministic  design  are 
shown  in  solid  lines,  while  the  dashed  lines  indicate  the  membership  functions  of  the  fuzzy  set 
based  design.  The  possibility  distributions  of  Fig.  (5.2)  illustrate  the  differences  in  the  manner 
each  approach  maximizes  a safety  measure  for  a given  weight,  as  discussed  in  Section  5.1.  For 
the  deterministic  optimum  design,  the  solid  lines  converge  to  a single  point  for  a possibility  of 
failure  equal  to  1,  indicating  that  the  deterministic  approach  attempts  to  equalize  the  factor  of 
safety  associated  with  the  Von  Mises  yield  stress  failure  criterion  and  that  associated  with  the 
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Figure  5.2:  Possibility  distribution  of  failure  for  the  deterministic  and  fuzzy  set  based 

optimum  designs. 

buckling  load  constraint.  For  the  fuzzy  set  based  optimum  design,  the  dashed  lines  fall  within 
the  solid  lines  at  (P  - P/)=0,  indicating  that  the  fuzzy  set  based  approach  attempts  to  equalize 
the  possibility  of  failure  values  due  to  yielding  and  due  to  violating  the  buckling  load 
constraint.  However,  for  the  fuzzy  set  based  optimum  design  the  design  variable  limits  kept 
the  two  possibility  of  failure  values  from  becoming  equal. 

The  factor  of  safety,  calculated  for  the  optimum  fuzzy  set  based  design,  is  not  much 
different  from  that  of  the  deterministic  design  (only  0.7%  lower),  however,  there  is  a large 
difference  in  the  possibility  of  failure  between  the  two  designs.  The  maximum  possibility  of 
failure  for  the  fuzzy  set  based  design  is  22.3%  lower  than  that  of  the  deterministic  design. 

5,4.  Cost  versus  Weight  Design  Chart 

Finally,  the  dependence  of  the  weight  of  the  structure  on  the  level  of  uncertainty 
associated  with  the  design  variables  was  also  investigated.  For  this  investigation,  the  allowable 
possibility  of  failure  was  kept  constant  at  the  maximum  possibility  of  failure  for  the  optimum 
fuzzy  set  based  design  obtained  from  Eqs.  (5.5)  and  (5.6)  (i.e.,  0.2788  as  shown  in  Table  5.3). 
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The  level  of  uncertainty  associated  with  the  design  variables  A,  p and  y,  was  varied  between 
± 2%  and  ± 20%.  As  for  the  fuzzy  set  based  design  problem,  the  reduced  possibility  of  failure 
response  surface  approximations  of  CHAPTER  4 were  used  to  calculate  the  possibility  of 
failure. 

5.4.1. Optimization  Problem 

For  constructing  the  weight  versus  cost  design  chart,  seven  evenly  distributed  levels  of 
the  uncertainty,  u,  between  ± 2%  and  ± 20%  were  considered.  For  each  of  these  levels  the 
non-dimensional  cross-sectional  area  of  the  plate  was  minimized  for  the  specified  allowable 
possibility  of  failure.  The  corresponding  optimization  problem  may  then  be  written  as: 
Minimize: 

^"^_"|(1  + 2/?  + r)  + f(1“2/?_r)  (5-7) 

Subject  to: 


0.4 


1>0 


0.4 


y > 0 


{ r + P 

0.4 


>o 


^(p-p /?,/,«) 

0.2788 


-1>0 


(5.8) 


5.4.2. Results 

Recall  that  it  is  assumed  that  the  manufacturing  cost  of  the  structure  is  dominated  by 
the  level  of  uncertainty  in  the  design  variables,  while  the  non-dimensional  cross-sectional  area 
is  representative  of  the  weight.  Increasing  the  level  of  uncertainty  in  the  design  variables 
corresponds  to  an  increase  in  the  allowable  tolerances  in  these  variables  and  would  thus 
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correspond  to  a reduction  in  manufacturing  cost.  In  contrast,  an  increase  in  the  non- 
dimensional  cross-sectional  area  would  correspond  to  an  increase  in  weight. 

As  expected,  the  non-dimensional  cross-sectional  area  of  the  plate  increased  with  an 
increase  in  the  level  of  uncertainty  and  the  results  are  shown  graphically  in  Fig.  (5.3). 
Figure  (5.3)  was  constructed  based  on  seven  optimizations,  indicated  by  the  circular  symbols. 
Figure  (5.3)  indicates  that  the  increase  in  weight  is  almost  proportional  to  the  increase  in  the 
uncertainty  associated  with  the  design  variables.  The  non-dimensional  cross-sectional  area 
increased  by  10.7%  with  an  18%  increase  in  the  uncertainty  associated  with  the  design 
variables.  Using  Fig.  (5.3)  and  the  dependence  of  the  manufacturing  cost  on  the  tolerance  of 
A,  j5  and  the  designer  may  determine  what  tolerance  to  use  when  manufacturing  the  plate, 
by  trading  of  manufacturing  cost  with  weight. 


Level  of  Uncertainty  in  Design  Variables  [xlOO  %] 


Figure  5.3:  Non-dimensional  cross-sectional  area  compared  to  the  level  of  uncertainty  in 

the  design  variables. 


5.5.  Reduction  in  Computational  Cost 

Response  surface  approximations  formed  an  integral  part  of  the  fuzzy  set  based  design 
process  presented  here.  The  reduced  stress  distribution  and  buckling  load  response  surface 
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approximations  of  CHAPTER  4 were  used  to  replace  computationally  expensive  finite 
element  analyses  during  the  evaluation  of  the  possibility  of  failure.  Fuzzy  set  based  designing 
for  uncertainty  is  an  expensive  two  level  optimization  problem  and  the  use  of  the  response 
surface  approximations  resulted  in  substantial  reductions  in  computational  costs  due  to  the 
following  two  reasons: 

1.  It  is  cheaper  to  construct  all  of  the  required  approximations  than  to  perform  a single 
two  level  optimization  problem. 

2.  The  computational  cost  is  shifted  from  the  optimization  problem  to  that  of 
constructing  the  approximations.  Shifting  the  computational  cost  to  the  problem  of 
constructing  the  approximations  allows  the  use  of  the  approximations  in  additional 
optimization  problems  at  very  little  additional  computational  cost. 

The  reductions  in  computational  cost  are  illustrated  for  a single  two  level  optimization 
problem,  by  considering  all  problem  parameters  to  be  uncertain,  as  is  the  case  in  the  present 
design  problem.  The  evaluation  of  the  possibility  of  failure  for  a single  a level  cut  then 
requires  2*28=512  finite  element  analyses  when  no  response  surface  approximations  are  used. 
For  a single  optimization,  a conservative  estimate  of  the  required  number  of  finite  element 
analyses  required,  when  not  using  response  surface  approximations,  is  obtained  from  the 
product  of  the  following  four  numbers: 


Average  number  of  design  optimization  iterations:  5 

Average  number  of  n(P.P/)  evaluations  per  iteration:  6 

Average  number  of  a level  cut  evaluations  per  n(P.P/)  evaluation:  5 

Number  of  finite  element  analyses  per  a level  cut  evaluation  of  11(P.P^):  512 

Total  number  of  finite  element  analyses  required  per  optimization: 


76,800 
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In  contrast,  the  stress  distribution  and  buckling  load  response  surface  approximations 
were  constructed,  and  evaluated,  from  a total  of  only  752  finite  element  analyses.  These 
results  indicate  that  response  surface  approximations  provide  an  effective  approach  for 
reducing  the  computational  cost  associated  with  performing  fuzzy  set  based  design.  The  large 
number  of  computationally  expensive  finite  element  analyses  required  to  perform  the  fuzzy 
set  based  design  is  replaced  by  response  surface  approximations  that  are  inexpensive  to 
evaluate. 

By  using  response  surface  approximations,  the  computational  burden  shifts  from  the 
optimization  problem  to  the  problem  of  constructing  the  response  surface  approximations. 
Due  to  the  iterative  nature  of  the  design  process,  the  fact  that  response  surface  approximations 
allow  multiple  optimizations  at  minimal  additional  computational  cost  should  be  an  attractive 
feature  to  any  designer.  Additionally,  the  fact  that  multiple  optimizations  may  be  performed 
at  little  additional  computational  cost  allows  the  designer  to  investigate,  for  example,  the 
dependence  of  the  weight  on  the  uncertainty  in  the  design  variables  in  the  form  of  weight 
versus  cost  design  charts.  The  reduced  possibility  of  failure  response  surface  approximation 
was  used  mainly  to  simplify  the  integration  of  the  analysis  code  and  the  optimization 
algorithm,  but  also  to  filter  out  numerical  noise  by  providing  an  approximate  response 
function,  thus  allowing  the  use  of  a derivative-based  optimization  algorithm. 


CHAPTER  6 

DROPPED  PLY  COMPOSITE  LAMINATED  PLATE  RESPONSE  SURFACE 

CONSTRUCTION 

Designing  for  uncertainty  of  a composite  laminated  plate  with  internally  terminated 
plies  (referred  to  as  dropped  plies)  is  the  second  example  problem  considered.  The  dropped 
plies  result  in  a change  in  thickness  across  the  width  of  the  plate  and  are  used  to  reduce  the 
weight  of  the  structure.  The  weight  is  reduced  by  reducing  the  thickness  according  to  the 
loads  so  that  the  structure  is  thicker  in  regions  with  higher  loads.  It  is  assumed  that  the 
structure  will  be  built  from  future  materials,  with  little  information  available  on  uncertainty 
in  the  material  properties,  which  will  be  modeled  by  fuzzy  set  theory.  The  present  problem  is 
thus  an  extension  of  the  isotropic  plate  example  problem,  where  the  complexity  of  the 
problem  is  increased  by  considering  a composite  laminated  plate.  For  the  present  problem  the 
material  properties  and  applied  load  are  considered  to  be  uncertain.  Additionally, 
manufacturing  tolerances  are  also  modeled  by  considering  the  uncertainty  in  the  lay-up  angle 
of  the  individual  plies  of  the  composite  laminated  plate.  The  design  optimization  trades  off 
costs  associated  with  material  property  characterization  versus  costs  associated  with  increased 
structural  weight. 

6.1.  Problem  Description 

Reducing  the  weight  by  reducing  the  thickness  of  composite  laminated  structures 
according  to  the  loads  is  common  in  the  aircraft  industry.  Examples  include  helicopter  rotor 
blades  as  well  as  wing  structures,  where  the  thickness  is  reduced  from  the  root  to  the  tip.  The 
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internal  termination  of  plies  is  an  effective  method  for  reducing  the  thickness  of  composite 
laminated  structures.  Unfortunately,  dropped  plies  result  in  severe  stress  concentrations  due 
to  the  stiffness  and  thickness  discontinuities  associated  with  the  internal  termination  of  the 
plies.  These  stress  concentrations  may  greatly  reduce  the  strength  of  the  structure,  thus 
offsetting  the  advantages  of  reducing  the  thickness  by  terminating  plies. 

The  dropped  ply  problem  has  been  studied  extensively,  and  a large  number  of 
analytical,  numerical  and  experimental  studies  has  been  published  since  the  mid  1980’s.  For 
example  Kemp  and  Johnson77  (1985)  performed  a numerical  study  on  graphite-epoxy 
composite  laminated  plates  subjected  to  both  a tensile  and  a compressive  applied  load,  while 
Fish  and  Lee80  (1988)  performed  a numerical  and  experimental  study  on  tensile  loaded  glass- 
epoxy  composite  laminated  plates.  Dinardo  and  Lagace81  (1989)  performed  both  an  analytical 
and  experimental  study  of  the  buckling  and  post-buckling  behavior  of  graphite-epoxy 
composite  laminated  plates,  while  Curry,  Johnson  and  Starnes82  (1992)  performed  a numerical 
and  experimental  study  of  similar  structures  subjected  to  both  a tensile  and  a compressive 
applied  load.  Additionally,  Llanos  and  Vizzini83  (1992)  investigated  the  influence  of  adding 
film  adhesive  on  the  strength  of  tension  loaded  glass-epoxy  specimens,  while  Fish  and 
Vizzini84,8'  (1992  and  1993)  Botting,  Vizzini  and  Lee86  (1996)  studied  the  influence  of  changing 
the  ply  drop  configuration.  Davila  and  Johnson87  (1993)  performed  a combination  of 
numerical  and  experimental  analyses  of  the  buckling  and  post-buckling  behavior  of  graphite- 
epoxy  composite  laminated  plates  and  Harrison  and  Johnson88  (1993)  and  Rose  and  Starnes89 
(1996)  formulated  efficient  approximate  analyses  for  determining  the  inter-laminar  stresses  in 
composite  laminated  plates  with  thickness  discontinuities.  Finally,  Vizzini90  (1995)  studied  the 
influence  of  including  imperfections  that  naturally  occur  during  the  manufacturing  process,  in 
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the  modeling  of  glass-epoxy  structures.  The  experimental  studies  indicated  that  delamination 
is  a common  mode  of  failure.  As  a result,  all  of  the  above  mentioned  studies  concentrated  on 
failure  initiation  in  the  form  of  delamination  that  occurs  in  the  interfaces  directly  above  and 
below  the  dropped  plies.  The  studies  indicated  that  the  presence  of  dropped  plies  greatly 
reduces  the  strength  of  the  resulting  composite  laminated  plate. 

In  the  present  example  problem,  the  cost  design  of  a composite  laminated  plate  with 
dropped  plies  is  considered.  The  plate  structure  is  typical  of  wing  skin  applications,  where 
aerodynamic  requirements  for  a smooth  surface  introduce  asymmetry.  This  problem  was 
studied  extensively  by  Curry,  Johnson  and  Starnes82  (1992)  and  by  Davila  and  Johnson87  (1993) 
and  is  shown  schematically  in  Fig  (6.1). 


Figure  6.1:  Dropped  ply  composite  laminated  plate  typical  of  wing  skin  applications. 

Curry,  Johnson  and  Starnes82  (1992)  performed  an  experimental  study,  considering  37 
specimens  that  differed  only  in  the  configuration  of  the  dropped  plies.  The  experimental 
study  revealed  that  the  axial  strength  of  the  dropped  ply  composite  laminated  plate  considered 
is  less  than  that  of  the  thin  section,  and  that  the  compressive  strength  is  less  than  the  tensile 
strength.  The  authors  also  concluded  that  the  reduction  in  axial  strength  is  directly  related  to 
the  change  in  axial  stiffness  between  the  thick  and  the  thin  sections. 
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Davila  and  Johnson87  (1993)  correlated  the  experimental  results  of  Ref.  (82)  with 
detailed  finite  element  analyses  and  proposed  a delamination  initiation  failure  index  based  on 
an  existing  quadratic  inter-laminar  failure  index.  The  proposed  failure  index  is  evaluated  at  a 
characteristic  distance  from  the  ply  drop  stress  singularity  and  correlated  well  with  the  failure 
initiation  observed  in  the  experimental  data  of  Ref.  (82).  However,  the  success  of  the  method 
depends  on  the  availability  of  experimental  data  for  determining  the  characteristic  distance. 

The  present  problem  is  based  on  the  254.0  mm  by  25.46  mm,  AS/3502  graphite-epoxy 
(material  properties  summarized  in  Table  6.1)  plate  shown  in  Fig.  (6.1),  that  was  considered  by 
Curry,  Johnson  and  Starnes8'  (1992)  and  by  Davila  and  Johnson87  (1993).  However,  in  the 
present  problem  the  plies  are  dropped  in  two  steps,  separated  by  20  ply  thicknesses,  and  the 


Table  6.1:  Composite  laminated  plate  material  properties. 


AS4/3502  Graphite-Epoxy 

E, 

126  GPa 

Young’s  modulus  in  fiber  direction 

e2 

11.3  GPa 

Young’s  modulus  in  matrix  direction 

E3 

11.3  GPa 

Young’s  modulus  in  thickness  direction 

Vn 

0.3 

Poisson’s  ratio  in  1-2  plane 

Vt3 

0.3 

Poisson’s  ratio  in  1-3  plane 

V '23 

0.35 

Poisson’s  ratio  in  2-3  plane 

gI2 

6.0  GPa 

Shear  modulus  in  1-2  plane 

g13 

6.0  GPa 

Shear  modulus  in  1-3  plane 

g23 

3.38  GPa 

Shear  modulus  in  2-3  plane 

Neat  Resin 

E 

3.45  GPa 

Young’s  modulus 

V 

0.41 

Poisson’s  ratio 

G 

1.34  GPa 

Shear  modulus 

Strength  Allowables 

xt 

1.45  GPa 

Fiber  direction  tensile  strength 

Xc 

-1.45  GPa 

Fiber  direction  compressive  strength 

Y„Z, 

52  MPa 

Transverse  tensile  strength 

Yc  Zc 

-206  MPa 

T ransverse  compressive  strength 

Si2,  SI3,  S23 

93.1  MPa 

Shear  strengths 
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resin  pockets  are  reinforced  with  90°  fibers,  measure  relative  to  the  length  direction  of  the 
plate  (i.e.,  the  x-direction).  The  present  problem  is  shown  schematically  in  Fig.  (6.2).  The 
drop  zone  geometry  used  in  the  present  problem  was  determined  based  on  an  investigation 
using  finite  element  analyses  as  discussed  in  Appendix  A. 


Normal  (Z) 


Transverse  (V) 
Longitudinal  (X) 


Continuous  Plies 
(Quasi-Isotropic) 

Drop  Plies 
(0°  Plies) 

Continuous  Plies 
(Quasi-Isotropic) 


Thick  Section  Transition  Section  Thin  Section 


First  Ply  Drop 


Second  Ply  Drop 


Figure  6.2:  Two  step  dropped  ply  composite  laminated  plate. 

As  shown  in  Fig.  (6.2),  the  laminated  plate  consists  of  three  sub-laminates  where  four 
0°  plies  are  dropped  and  the  sub-laminates  surrounding  the  dropped  plies  are  quasi-isotropic. 
More  specifically,  the  lay-up  of  the  three  sub-laminates  in  the  thick  section  of  the  plate  is 
[(±45/0/90)s,(02)s,(  ±45/0/90)s]t.  The  three  sub-laminates  are  reduced  to  two  sub-laminates  in 
the  thin  section,  with  a lay-up  of  [(±45/0/90)s,(  ±45/0/90)s]T.  A detailed  cross-section  of  the 
dropped  ply  composite  laminated  plate  is  shown  in  Fig.  (6.3). 

Two  failure  mechanisms  were  considered  in  the  design,  delamination  in  the  interfaces 
directly  above  and  below  the  dropped  plies,  and  ply  failure.  Response  surface  approximations 
of  the  stress  and  strain  components  required  to  evaluate  these  failure  mechanisms  (discussed 
subsequently  in  Section  6.2)  were  constructed  and  used  during  the  calculation  of  the 
possibility  of  failure  of  the  composite  laminated  plate.  These  response  surface  approximations 
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Normal  (Z) 

* 


Quasi-Isotropic 
(8  Plies) 

All  0°  Plies 
(4  Plies) 

Quasi-Isotropic 
(8  Plies) 


0.5028  mm  0.5028  mm 


Longitudinal  (X)  (4tpty) 

H H 


127.000  mm 

I 2.514  mm 

124.486  mm 

1 (20  tD.v) 

Figure  6.3:  Cross-section  model  with  dimensions  shown. 

were  constructed  from  numerical  experiments  in  the  form  of  linear  finite  element  analyses 
with  CSAR/NASTRAN  Version  9 791  (1997). 

Due  to  symmetry,  only  half  the  width  (i.e.,  y-direction)  of  the  plate  was  modeled  and 
the  finite  element  model  had  a total  of  2,814  elements  and  4,020  nodes.  The  model  consisted 
of  a combination  of  four-node  plate  bending  elements  and  eight-node  solid  elements  as  shown 
schematically  in  Fig.  (6.4).  Multi-point  constraints  were  used  to  connect  the  plate  bending 
elements  (with  rotational  degrees  of  freedom)  to  the  solid  elements  (without  rotational  degrees 
of  freedom). 


(a)  (b) 

Figure  6.4:  Schematic  finite  element  model  (actual  model  had  42  two-dimensional  and 

2,772  three-dimensional  elements):  (a)  Overall  perspective;  (b)  Detailed 

perspective  near  the  ply  drop. 
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The  four-node  plate  bending  elements  were  used  to  model  the  plate  away  from  the  ply 
drop,  while  the  eight-node  solid  elements  were  used  for  the  detailed  modeling  of  a small  region 
near  the  ply  drop.  Three  dimensional  elements  are  required  in  the  vicinity  of  the  ply  drop  in 
order  to  model  the  three  dimensional  stress  field  in  this  region.  Although  multiple  three 
dimensional  elements  were  used  through  the  thickness  of  each  sub-laminate,  these  elements  do 
not  correspond  to  the  ply  lay-up  of  the  sub-laminates.  Instead  each  sub-laminate  was  defined 
in  terms  of  its  equivalent  homogeneous  orthotropic  material  properties.  The  equivalent 
homogeneous  orthotropic  material  properties  were  calculated  using  the  methodology  as 
outlined  in  Sun  and  Li92  (1988). 

The  boundary  conditions  of  the  finite  element  model  were  chosen  to  approximate  the 
experimental  setup  used  by  Curry,  Johnson  and  Starnes82  (1992)  and  are  shown  schematically 
in  Fig.  (6.5).  A compressive  load  of  507.3  N/mm  was  applied  at  the  thick  end  (i.e.,  x=0  mm) 
where  only  a uniform  displacement  in  the  x-direction  was  allowed  (all  other  degrees  of 
freedom  were  fixed).  At  the  thin  end  (i.e.,  x=254.0  mm)  all  six  degrees  of  freedom  were  fixed. 


Reaction  Edge  Conditions 
U = V = W = 5W/5X  = 0 


Loaded  Edge  Conditions 
U = Constant 
V = W=5W/5X=0 


Figure  6.5:  Schematic  representation  of  the  boundary  conditions. 
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Additionally,  the  out  of  plane  displacements  were  set  equal  to  zero  within  a 9.525  mm  region 
along  the  width  of  each  endpoint  (i.e.,  x=0  mm  and  x=254.0  mm).  Along  the  length  of  the 
plate,  the  structure  was  supported  by  two  knife  edges  applied  3.175  mm  inward  from  the  free 
edges  (i.e.,  y=0  mm  and  y=25.46  mm).  For  the  shell  elements,  the  knife  edges  were  modeled 
by  setting  the  out  of  plane  displacement  equal  to  zero.  For  the  solid  elements,  the  knife-edges 
were  modeled  by  setting  the  out  of  plane  displacement  equal  to  zero  along  the  bottom  of  the 
plate  (i.e.,  z= 0 mm).  Additionally,  a small  pinching  displacement  equal  to  -lxlO-4  tQ  , where  ta 
is  the  thickness  of  the  thick  side  of  the  plate,  was  imposed  along  the  top  surface  of  the  thick 
section  of  the  plate.  This  pinching  displacement  was  used  to  simulate  the  clamping  action  of 
the  knife  edge  support. 

6.2.  Analysis  Issues 

The  baseline  structure  (nominal  values  of  problem  parameters  as  shown  in  Table  6.1) 
was  investigated  in  order  to  identify  the  critical  region  for  each  of  the  delamination  and  ply 
failure  mechanisms.  Response  surface  approximations  of  the  stress  and  strain  components 
required  to  evaluate  the  failure  mechanisms  were  then  constructed  in  the  critical  regions 
identified  by  this  investigation. 

Based  on  work  by  Whitney  and  Nuismer93  (1974),  a characteristic  distance  was  used  to 
account  for  the  stress  singularities,  and  both  the  delamination  and  ply  failure  criteria  were 
evaluated  at  a characteristic  distance  from  the  corresponding  stress  concentration.  The 
characteristic  distance  must  be  determined  experimentally.  The  results  of  Davila  and 
Johnson87  (1993)  indicated  that  the  characteristic  distance  depends  on  both  the  material  used 
and  the  lay-up  of  the  composite  laminated  plate.  Due  to  a lack  of  experimental  data  in  the 
present  dissertation,  it  was  assumed  that  both  failure  criteria  are  evaluated  at  an  equal 
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characteristic  distance  of  0.25  mm.  This  value  of  the  characteristic  distance  was  based  on 
results  obtained  by  Davila  and  Johnson87  (1993)  for  the  problem  shown  in  Fig.  (6.1).  Note 
that  neither  of  the  assumptions  (i.e.,  equal  characteristic  distance  for  the  two  failure  modes  and 
a characteristic  distance  of  0.25  mm)  restricts  the  generality  of  the  proposed  methodology. 
The  response  surface  approximations  of  the  stress  and  strain  distributions,  used  to  calculate 
the  possibility  of  failure,  do  not  depend  on  the  characteristic  distance. 

6.2.1. Ply  Failure 

The  maximum  stress  failure  criterion  (Gibson94  1994,  pp.  103-106)  was  used  to  find  the 
most  critical  region  for  ply  failure.  The  maximum  stress  criterion  predicts  failure  when  any 
stress  component  in  the  principal  material  system  exceeds  the  corresponding  strength,  and 
thus  to  avoid  failure  the  following  inequalities  must  be  satisfied 

< cr,  < X,  Yc  < cr2  < Y,  |cr12 1 < Su  (6.1) 

where  the  subscript  1 indicates  the  longitudinal  fiber  direction  and  the  subscript  2 indicates 
the  transverse  fiber  direction  as  shown  in  Fig.  (6.6).  Additionally,  the  subscript  c denotes 
compression  and  the  subscript  t denotes  tension. 


Figure  6.6:  Composite  laminated  plate  fiber  directions. 
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The  maximum  stress  criterion  must  be  applied  to  each  ply  since,  unlike  homogenous, 
isotropic  materials,  it  is  not  obvious  which  is  the  critical  ply.  In  the  present  dissertation,  the 
topmost  ply  on  the  thinnest  section  of  the  laminated  plate  was  identified  as  the  critical  ply. 
This  ply  corresponds  to  the  -45°  ply  ahead  of  the  second  ply  drop  as  shown  in  Fig.  (6.2).  For 
this  ply  cr,  and  cr2  are  in  compression  and  the  failure  indices,  based  on  the  maximum  stress 
criterion  of  Eq.  (6.1)  are  defined  as: 


(6.2) 


Failure  thus  occurs  if  any  of  the  failure  indices  of  Eq.  (6.2)  is  equal  to  or  larger  than  1,  with 
higher  failure  index  values  indicating  more  critical  stress  components.  The  failure  indices 
calculated  from  the  baseline  problem,  with  no  uncertainty  in  the  problem  parameters,  are 
shown  graphically  in  Fig.  (6.7)  for  the  topmost  ply  of  each  of  the  two  thinner  sections  of  the 
composite  laminated  plate.  These  thinner  sections  correspond  to  the  first  and  the  second  ply 
drops  as  shown  in  Fig.  (6.2). 


Distance  Measured  From  First  Ply  Drop  [mm] 


Figure  6.7:  Ply  failure  indices  for  the  topmost  ply  of  the  two  thinner  sections. 
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In  Fig.  (6.7)  the  x-axis  represents  a local  coordinate  system  with  origin  at  the  first  ply 
drop  and  the  solid  vertical  lines  indicate  the  characteristic  distance.  Figure  (6.7)  indicates  that 
the  shear  stress  component  on  the  top  surface  of  the  thinnest  section  of  the  dropped  ply 
composite  laminated  plate  is  the  most  critical  stress  component.  To  be  more  consistent  with 
the  convention  used  by  industry,  the  maximum  strain  failure  criterion  was  used  throughout 
the  design  process,  by  converting  the  shear  stress  allowable  to  a strain  allowable  given  by: 


(6.3) 


The  critical  ply  failure  index  FPiy  of  the  structure  may  then  be  written  in  terms  of  the 
Yn  shear  strain  on  the  top  surface  of  the  thinnest  section  of  the  laminated  plate  and  the 
allowable  shear  strain  as  follows: 


p -hL 

1 Ply  ~ 

en 


(6.4) 


The  ply  failure  load  PPIy  is  the  defined  as 


Ppiy 


Ply 


(6.5) 


where  P denotes  the  applied  load.  This  failure  index  is  evaluated  at  the  characteristic  distance 
away  from  the  re-entrant  corner. 


6.2.2. 


Delamination  Failure 


For  the  delamination  failure,  the  failure  index  proposed  and  successfully  evaluated  by 
Davila  and  Johnson87  (1993)  was  used.  This  failure  index  is  based  on  a quadratic  failure  index 
introduced  by  Brewer  and  Lagace95  (1988) 
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where  crzy,  aZY  and  crzz  are  the  inter-laminar  stress  components  and  SI2,  S2 3 and  Zt  are  the  inter- 
laminar strength  allowables  (see  Table  6.1).  This  failure  index  was  evaluated  in  the  interfaces 
directly  above  and  below  the  dropped  plies,  at  a characteristic  distance  away  from  the  ply 
drop.  Additionally,  in  order  to  avoid  free  edge  failure,  the  failure  index  was  only  evaluated 
inwards  from  the  knife  edge  support  shown  in  Fig.  (6.5).  The  delamination  failure  load  PDeiam 
is  calculated  from  the  failure  index  as: 


Delam 


(6.7) 


To  identify  the  critical  region  for  delamination  failure,  the  delamination  failure  index 
pDeiam  was  evaluated  for  the  baseline  problem  with  no  uncertainty  in  the  problem  parameters 
and  the  results  are  shown  graphically  in  Fig.  (6.8).  In  Fig.  (6.8)  the  x-axis  represents  a local 
coordinate  system  with  origin  at  the  first  ply  drop  and  the  solid  vertical  lines  indicate  the 
characteristic  distance.  Figure  (6.8)  indicates  that  the  interface  above  the  second  step  of  the 
ply  drop  corresponding  to  the  transition  section  as  shown  in  Fig.  6.2  is  the  most  critical. 


Distance  Measured  From  First  Ply  Drop  [mm] 


Figure  6.8: 


Delamination  failure  index  above  and  below  the  dropped  plies. 
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6.3.  Dimensional  Analysis 

As  for  the  isotropic  plate  example  problem,  dimensional  analysis  was  used  to  reduce 
the  number  of  predictor  variables  of  the  stress  and  strain  distribution  response  surface 
approximations,  where  the  predictor  variables  are  the  variables  used  in  the  response  surface 
approximations.  For  the  present  example  problem,  the  material  properties,  characteristic 
distance,  lay-up  angles  and  applied  load  were  considered  as  uncertain  and  the  stress  and  strain 
distributions  were  thus  required  as  a function  of  these  uncertain  problem  parameters.  Unlike 
the  isotropic  plate  example  problem,  the  predictor  variables  of  the  present  problem  were  thus 
known  beforehand  and  the  Buckingham  n Theorem  (see  Appendix  B)  could  be  used  to  non- 
dimensionalize  the  problem,  using  the  functional  form  of  the  stress  and  strain  distributions. 

The  stress  and  strain  equations  may  be  written  in  functional  form  in  terms  of  the 
material  properties,  the  geometry  and  the  boundary  conditions  and  applied  load.  The 
fundamental  dimensional  units  of  both  the  stress  and  strain  distributions  are  force  [F]  and 
length  [L].  Thus,  both  the  stress  and  strain  distributions  may  be  written  in  non-dimensional 
form  by  selecting  E,  (dimension  [F.L'2])  and  a dummy  variable  / of  unit  length  (dimension  [L]) 
as  the  predictor  variables  with  independent  dimensions. 

Since  a linear  finite  element  analysis  is  performed,  both  the  stress  and  the  strain 
distributions  depend  linearly  on  the  magnitude  of  the  applied  load.  The  applied  load  can  thus 
be  eliminated  as  a predictor  variable,  since  the  stress  and  strain  distributions  can  simply  be 
scaled  for  different  values  of  the  applied  load  as  follows: 

P P 

r = Yo—  <y  = <7  — 


(6.8) 
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In  Eq.  (6.8)  P0  is  the  applied  load  used  in  the  liner  finite  element  analysis  and  cr0  and  ya  are  the 
corresponding  stress  and  strain  values,  while  P represents  a different  value  of  the  applied  load 
and  a and  y are  the  corresponding  stress  and  strain  values. 

Additionally,  because  equivalent  orthotropic  material  properties  were  used  in  the 
modeling  of  each  sub-laminate,  the  number  of  variables  required  to  characterize  the  lay-up  of 
each  sub-laminate  may  be  greatly  reduced  by  making  use  of  lamination  parameters.  Using 
lamination  parameters,  the  lay-up  of  each  sub-laminate  may  be  characterized  by  only  two  non- 
dimensional  lamination  parameters  (Haftka  and  Giirdal18  pp.  415-429)  instead  of  the  lay-up 
angle  of  each  ply,  as  follows 

v, * = -fJCos(20l)  v;  = -tcos(  40,)  (6.9) 

n i=o  n ,=o 


where  it  is  assumed  that  all  plies  have  the  same  thickness.  Additionally,  n denotes  the  total 
number  of  plies  in  the  sub-laminate  and  0 the  lay-up  angle  of  the  ith  ply. 

Based  on  the  above  observations,  the  non-dimensional  functional  relationships  of  the 
stress  and  strain  distributions  may  be  written  in  terms  of  the  predictor  variables  as 
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(6.10) 


where  the  subscripts  Drop,  Top  and  Bottom  are  used  to  denote  the  respective  sub-laminates. 
Additionally,  it  was  assumed  that  E^E},  Gi2=G3,  and  Oij=Vi3  as  shown  in  Table  (6.1). 
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Using  the  non-dimensional  lamination  parameters,  the  number  of  predictor  variables 
required  to  characterize  the  lay-up  of  each  sub-laminate  is  greatly  reduced  from  the  lay-up 
angle  of  each  ply  (total  of  20  variables)  to  only  two  lamination  parameters  per  sub-laminate 
(total  of  6 variables).  Additionally,  the  remaining  predictor  variables  were  reduced  from  14  to 
12  by  making  use  of  dimensional  analysis  and  the  fact  that  the  stress  and  strain  distributions 
depend  linearly  on  the  applied  load. 

6.4.  Calculating  the  Possibility  of  Failure 

The  possibility  ol  failure  is  calculated  from  the  load  margin  (P  - P/)  of  the  structure, 
with  failure  defined  to  occur  when 

P-Pf>  0 (6.11) 

and  the  failure  load  Pf  is  calculated  from  the  ply  and  delamination  failure  loads  as  follows: 

Pf  (6.12) 

b Delam 

For  the  present  example  problem,  the  design  space  is  not  a rectangular  domain  in  all  of 
the  design  variables.  In  particular,  the  design  space  for  the  lamination  parameters  V‘  and  V’ 
of  the  top  and  bottom  sub-laminates  are  non-rectangular  and  are  shown  graphically  in 
Fig.  (6.9)  for  a 5 degree  uncertainty  in  the  ply  lay-up  angle. 

Note  that  there  exists  a linear  relationship  between  V'  Drop  andF3*flrop.  The  total 

number  of  predictor  variables  is  thus  reduced  from  12  to  only  11  by  eliminating  F3*  Drop  as  a 
predictor  variable.  Apart  from  the  V*  and  F3*  design  spaces  of  the  top  and  bottom  sub- 
laminates that  are  non-rectangular,  both  the  V'  Drop  design  space  and  the  V'  and  F,"  design 
spaces  for  the  top  and  bottom  sub-laminates  are  also  discrete. 
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Figure  6.9:  All  possible  combinations  of  the  lamination  parameters  V’  and  V’ : 

(a)  Dropped  ply  sub-laminate;  (b)  Bottom  and  top  sub-laminates. 


The  Design  Optimization  Tools  (DOT)70  (1995)  software  was  used  to  verify  that  the 
global  optimum  of  the  load  margin  (P  - Pf)  (see  Eq.  6.11)  occurred  at  a vertex  of  the  design 
space,  thus  validating  the  use  of  the  vertex  method  to  calculate  the  possibility  of  failure.  Since 
the  sequential  linear  programming  algorithm  provided  with  the  Design  Optimization  Tools 
(DOT)70  (1995)  software  converges  on  a local  minimum  in  the  design  space,  1,000 
optimizations  were  performed  from  different  starting  points  randomly  distributed  throughout 
the  design  space  to  ensure  that  a global  optimum  is  found.  These  optimizations  were 
performed  for  the  zero  a level  cut  using  the  maximum  level  of  uncertainty  associated  with 
each  design  variable  (shown  subsequently  in  Table  6.3).  The  best  solutions  obtained  from  the 
1,000  optimizations  are  summarized  in  Table  (6.2).  For  the  delamination  load  margin  of 
Eq.  (6.11),  the  overall  minimum  value  was  found  on  the  boundary  of  the  V*  and  F3*  design 
spaces  for  the  top  and  bottom  sub-laminates,  but  not  at  a vertex.  Since  the  V’  and  V‘  design 
spaces  for  the  top  and  bottom  sub-laminates  are  discrete  (see  Fig.  6.9),  all  20  boundary  points 
for  each  of  the  two  V’  and  V'  design  spaces  were  included  in  the  vertex  method  to  ensure 
that  a global  optimum  is  found.  The  results  obtained  from  the  vertex  method  using  only  the 
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Table  6.2:  Maximal  load  margins  for  the  zero  a level  cut. 


Optimization 

Method 

P - P Ply 

P Poeiam 

Design  Optimization 
Tools  Software 

-4,130.79  N 

-6,351.83  N 

Vertex  Method 

-4,164.83  N 

(V ertices  Only) 
Vertex  Method 

-6,126.35  N 

-4,164.83  N 

(All  Boundary  Points) 

-6,365.14  N 

seven  vertex  points  for  each  of  the  Vx  and  V’  design  spaces  for  the  top  and  bottom  sub- 
laminates, as  well  as  the  results  obtained  from  the  vertex  method  using  all  20  boundary  points 
are  also  summarized  in  Table  (6.2). 

Note  that  the  sequential  linear  programming  algorithm  provided  with  the  Design 
Optimization  Tools  (DOT)70  (1995)  software  is  a continuous  optimization  algorithm,  while 
the  vertex  method  considers  only  a discrete  number  of  design  points.  However,  the  vertex 
method  gave  slightly  better  results  than  the  sequential  linear  programming  algorithm  (0.82% 
better  for  P - PPiy  and  0.21%  for  P - Poeiam  ) since,  as  expected,  the  sequential  linear 
programming  algorithm  converged  on  local  minima  in  the  design  space.  The  optimum  results 
for  (P 

P Delam  ) obtained  by  the  sequential  linear  programming  algorithm  are  shown 
graphically  in  Fig.  (6.10)  for  the  first  100  of  the  original  1,000  optimizations.  Ten  of  the  100 
optimizations  did  not  satisfy  all  of  the  constraints  within  the  allowable  number  of  design 
iterations,  which  was  set  at  500  iterations  due  to  computational  constraints.  Thus,  Fig.  (6.10) 
only  shows  the  optimum  results  of  the  remaining  90  optimizations.  The  results  shown  in 
Fig.  (6.10)  have  a 26.72%  variation  between  the  minimum  and  maximum  ( P - Poeiam  ) 
optimum  values,  indicating  the  need  for  an  expensive  global  optimization  algorithm,  like  the 
vertex  method,  to  calculate  the  possibility  of  failure  accurately. 
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Figure  6.10:  Optimal  results  obtained  by  the  Design  Optimization  Tools  software  for 

(P-  P Delam  ) • 


6.5.  Response  Surface  Approximation  Construction 
As  for  the  isotropic  plate  example  problem,  response  surface  approximations  were 
used  mainly  to  reduce  the  computational  burden  associated  with  calculating  the  possibility  of 
failure.  For  the  dropped  ply  composite  laminated  plate  problem,  the  failure  load  Pf  is 
calculated  from  (as  shown  in  Section  6.4) 

p, 

Delam 

and  failure  is  defined  to  occur  when: 

P-Pf>  0 (6.14) 

It  was  further  assumed  that  delamination  represents  catastrophic  failure  corresponding 
to  the  ultimate  load,  while  ply  failure  only  represents  the  onset  of  failure  corresponding  to  the 
limit  load.  In  the  aircraft  industry  the  ultimate  load  is  typically  assumed  to  be  1.5  times 
greater  than  the  limit  load,  which  is  incorporated  as  a factor  of  safety  for  the  delamination 
load,  while  a factor  of  safety  equal  to  1.0  was  used  for  the  ply  failure  load. 


(6.13) 
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As  shown  in  the  flow  diagram  of  Fig.  (4.5),  two  levels  of  response  surface 
approximations  were  constructed.  On  the  first  level,  accurate  response  surface 
approximations  were  constructed  for  evaluating  the  ply  failure  and  delamination  failure  loads. 
These  approximations  were  then  used  instead  of  computationally  expensive  finite  element 
analyses  during  the  calculation  of  the  possibility  of  failure  based  on  Eqs.  (6.13)  and  (6.14).  On 
the  second  level,  response  surface  approximations  were  constructed  for  the  possibility  of 
failure  as  a function  of  the  level  of  uncertainty  associated  with  the  uncertain  problem 
parameters.  This  approximation  was  constructed  mainly  to  simplify  the  integration  of  the 
analysis  code  and  the  optimization  algorithm.  However,  due  to  the  relatively  large  number  of 
uncertain  problem  parameters  in  the  present  example  problem,  calculating  the  possibility  of 
failure  from  the  delamination  and  ply  failure  response  surface  approximations  is 
computationally  intensive  in  itself.  The  second  level  response  surface  approximation  thus  also 
results  in  some  reduction  in  computational  cost. 

In  constructing  the  second  level  response  surface  approximation  for  the  possibility  of 
failure,  the  predictor  variables  were  obtained  by  dividing  the  uncertain  problem  parameters 
into  five  groups,  where  all  the  parameters  in  a specific  group  have  the  same  level  of 
uncertainty.  The  five  groups  are  the  elastic  material  properties,  the  ply  and  delamination 
failure  strength  allowables,  the  ply  lay-up  angle,  the  characteristic  distance  and  the  applied 
load.  The  predictor  variables  were  defined  as  the  uncertainty  associated  with  the  elastic 
material  properties  (EL),  the  ply  failure  strength  allowable  (SP),  the  ply  lay-up  angle  (PL)  and 
the  characteristic  distance  (CD).  The  uncertainty  associated  with  the  delamination  strength 
allowables  (SD)  was  also  considered,  but  not  as  an  independent  variable.  Instead,  the 
uncertainty  in  the  delamination  strength  allowables  was  assumed  to  be  twice  that  associated 
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with  the  ply  failure  strength  allowable.  Additionally,  the  weight  W of  the  structure  was 
considered  as  a fifth  predictor  variable  in  the  possibility  of  failure  response  surface 
approximation. 

The  uncertain  problem  parameters  and  the  nominal  level  of  uncertainty  associated 
with  each  parameter  are  summarized  in  Table  (6.3). 


Table  6.3:  Uncertain  problem  parameters. 


Property 

Description 

Uncertainty 

Group 

E, 

E2,  E3 
G12,  G/3 
G23 

Ol2,  O13 
023 

Young’s  Modulus 
Young’s  Modulus 
Shear  Modulus 
Shear  Modulus 
Poison’s  Ratio 
Poison’s  Ratio 

10% 

10% 

10% 

10% 

10% 

10% 

Elastic  Material 
Properties 
(EL) 

ei2 

Allowable  transverse  shear  strain 

20% 

Ply  Failure 

S,2 

Allowable  transverse  shear  strength 

40% 

Strength  Allowable 

Z , 

Allowable  transverse  tensile  strength 

40% 

(SP) 

a 

Lay-up  Anglen 

5° 

Ply  Lay-Up 

(PL) 

p 

Applied  load 

20% 

Load 

Characteristic 

d0 

Characteristic  Distance 

20% 

Distance 

m 

t The  uncertainty  associated  with  the  delamination  strength  allowables  are  assumed  to  be 
twice  as  large  as  that  associated  with  the  ply  failure  strength  allowable, 
ft  The  uncertainty  in  the  lay-up  angle  and  the  original  lay-up  of  each  sub-laminate  is  used  to 
calculate  Vt  and  V’  for  the  dropped,  bottom  and  top  sub-laminates  respectively. 


6.5.1. Stress  and  Strain  Distribution  Response  Surface  Approximations 

The  stress  and  strain  distribution  response  surface  approximations  have  a relatively 
large  number  of  predictor  variables  (11  predictor  variables)  and  the  computational  cost 
associated  with  evaluating  a single  finite  element  analysis  is  high.  Therefore,  statistical  design 
of  experiments,  in  particular  the  25-Optimality  criterion,  was  used  to  find  a small  set  of  data 
points  for  constructing  the  stress  and  strain  distribution  response  surface  approximations 
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without  compromising  accuracy.  The  22-Optimality  criterion  was  selected  since  the  design 
space  is  not  a rectangular  or  sphere-like  shape  in  all  the  design  variables  (see  Fig.  6.9),  and  the 
standard  designs  (e.g.,  Box-Behnken  and  Central  Composite  designs)  are  thus  not  applicable. 

A full  quadratic  response  surface  approximation  with  78  parameters  was  assumed  as  an 
initial  response  surface  approximation  for  the  stress  and  strain  distributions.  Note  that  each 
laminate  configuration  (corresponding  to  an  individual  finite  element  model)  is  described  by 
only  10  of  the  11  predictor  variables.  Each  finite  element  model  includes  a number  of  data 
points  with  different  x/l  values  (the  11th  predictor  variable)  corresponding  to  different  finite 
elements.  The  22-Optimality  design  criterion  was  used  to  select  200  data  points, 
corresponding  to  200  finite  element  analyses,  for  the  10  predictor  variables  that  influence  the 
configuration  of  the  composite  laminated  plate.  For  the  full  quadratic  response  surface 
approximation  with  78  parameters,  the  number  of  data  points  were  thus  chosen  to  be  roughly 
equal  to  two  and  a half  times  the  number  of  parameters  in  the  response  surface 
approximation.  The  number  of  data  points  identified  by  the  29-Optimality  criterion  is 
typically  selected  to  be  between  one  and  a half  and  three  times  the  number  of  parameters  in 
the  response  surface  approximation.  The  ratio  of  the  number  of  data  points  to  the  number  of 
parameters  is  heavily  influenced  by  the  cost  of  evaluating  the  response  at  a single  data  point. 
The  22-Optimal  data  points  were  selected  from  a set  of  57,049  candidate  points  that  were 
obtained  by  considering  a three  level  factorial  design  in  the  10  relevant  predictor  variables, 
consisting  of  three  evenly  spaced  levels  for  each  predictor  variable.  An  additional  50  data 
points,  randomly  distributed  throughout  the  design  space,  were  chosen  for  evaluating  the 
predictive  capabilities  of  the  response  surface  approximations. 
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For  each  finite  element  analysis,  the  stress  and  strain  values  at  all  node  points  within 
1 mm  of  the  relevant  stress  concentration  (i.e.,  the  ply  drop  for  delamination  failure  and  the 
re-entrant  corner  for  ply  failure)  were  recorded.  These  stress  and  strain  values  were  used  in 
constructing  the  response  surface  approximations  in  terms  of  the  11th  predictor  variable,  x/l. 
In  constructing  the  response  surface  approximations,  the  node  points  that  coincide  with  the 
stress  singularities  (i.e.,  the  node  points  at  the  ply  drop  and  the  re-entrant  corner)  were  ignored 
in  order  to  avoid  numerical  inaccuracies  associated  with  the  finite  element  models  at  these 
node  points.  The  finite  element  models  are  inaccurate  at  the  stress  singularities  since  the  stress 
and  strain  values  tend  to  infinity,  which  cannot  be  represented  accurately  by  a finite  element 
analysis.  Ignoring  the  node  points  at  the  stress  singularities  implies  that  the  response  surface 
approximations  would  only  predict  realistic  stress  and  strain  values  some  distance  away  from 
the  stress  singularities,  and  not  at  the  singularities  themselves.  Since  the  characteristic  distance 
used  in  evaluating  both  failure  criteria  spans  a number  of  node  points  in  the  finite  element 
model,  the  response  surface  approximations  were  sufficient  to  obtain  accurate  stress  and  strain 
values.  This  process  resulted  in  1,000  data  points  for  constructing  the  yn  strain  distribution 
and  1,200  data  points  for  constructing  each  of  the  Ozz,  <Jzx  and  aZy  stress  distribution  response 
surface  approximations  (Data  Set  1).  Additionally,  250  independent  data  points  were  available 
for  evaluating  the  predictive  capabilities  for  the  yn  strain  distribution,  and  300  independent 
data  points  were  available  for  evaluating  the  predictive  capabilities  for  each  of  the  Ozz,  <Tzr  and 
Gzy  stress  distribution  response  surface  approximations  (Data  Set  2). 

The  initial  quadratic  response  surface  approximation  was  used  to  identify  the  most 
important  variables  for  describing  the  response.  Each  variable  was  varied  between  its  upper 
and  lower  bounds  with  all  other  variables  kept  constant  and  the  resulting  influence  on  the 
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response  function  was  investigated.  The  variables  with  the  largest  influence  on  the  response 
function  were  identified  as  the  most  important  variables.  This  investigation  was  performed 
graphically  using  JMP67  (1995),  which  allows  the  user  to  easily  study  the  influence  on  the 
response  when  changing  a single  variable  for  different  levels  of  the  remaining  variables.  Four 
variables  were  identified  as  important  variables  and  all  cubic  parameters  in  these  variables 
were  added  to  the  initial  quadratic  response  surface  approximation.  Adding  these  cubic  terms 
increased  the  number  of  parameters  in  the  initial  response  surface  approximation  from  78 
to  98.  Note  that  a full  cubic  polynomial  would  result  in  a response  surface  approximation 
with  364  parameters.  Due  to  the  relatively  high  computational  cost  associated  with  evaluating 
the  response  at  a single  data  point,  evaluating  the  response  at  the  large  number  of  data  points 
required  for  constructing  a response  surface  approximation  with  364  parameters  was 
impractical. 

The  .D-Optimal  data  points  were  selected  based  on  the  original  quadratic  response 
surface  approximation  and  should  thus  be  used  with  caution  in  constructing  the  enhanced 
response  surface  approximation  that  includes  both  quadratic  and  cubic  parameters.  For  the 
present  example  problem,  the  number  of  iTOptimal  data  points  is  roughly  equal  to  two  times 
the  total  number  of  parameters  in  the  enhanced  model.  Additionally,  the  random  points  used 
for  evaluating  the  predictive  capabilities  of  the  response  surface  approximations  would 
indicate  when  the  selected  /^-Optimal  data  points  are  not  appropriate  for  constructing  the 
enhanced  response  surface  approximation.  When  the  selected  data  points  are  not  appropriate 
for  constructing  the  enhanced  response  surface  approximation,  the  error  terms  associated  with 
the  random  data  points  would  be  much  greater  than  the  error  terms  associated  with  the  data 
point  used  to  construct  the  response  surface  approximation. 
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6.5.1. 1.  Strain  distribution  response  surface  approximation  (Ply  failure) 

The  ply  failure  load  may  be  calculated  from  Eqs.  (6.4)  and  (6.5),  and  depends  only  on 
the  yu  strain  component  on  the  top  surface  of  the  thinnest  section  of  the  composite  laminated 
plate.  The  yI2  strain  component  may  be  written  in  non-dimensional  functional  form 
(see  Eq.  6.10)  as: 


Yn  = Yn 


r G^_  Gj3_  y y'  v-  y y 

F ^ F 9 F 9 12  ’ 23  ? X-Dr°P'>  X-T°P’>  '*-ToP’>  l _ Bottom’ r 3_Bottom’  , 

V i 


(6.15) 


Where  V3 ' _Drop  depends  linearly  on  V’Drop  and  is  not  required  as  a predictor  variable 
(see  Section  6.5) 

For  the  y2  strain  component,  the  quadratic  response  surface  approximation  was 
enhanced  by  including  all  third  order  parameters  for  the  four  most  important  variables,  which 
were  identified  as:  E2/Eh  G12/Eh  Vl  Drop  and  x/l.  The  number  of  parameters  in  the  enhanced 

response  surface  approximation  was  reduced  from  98  to  60  using  the  stepwise  regression 
procedure  and  the  Cp  statistic.  Using  both  Data  Set  1 and  Data  Set  2,  the  predictive 
capabilities  of  the  full  quadratic  response  surface  approximation  (78  parameters),  the  full 
enhanced  response  surface  approximation  (98  parameters),  and  the  reduced  enhanced  response 
surface  approximation  (60  parameters)  were  compared  and  the  results  are  summarized  in 
Table  (6.4). 

Table  (6.4)  indicates  that  the  enhanced  response  surface  approximations  are  more 
accurate  than  the  original  quadratic  response  surface  approximation,  with  the  reduced 
enhanced  response  surface  approximation  slightly  more  accurate  (except  for  the  %RMSE 
value)  than  the  full  enhanced  response  surface  approximation.  As  for  the  isotropic  plate 
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Table  6.4:  Strain  distribution  response  surface  approximation. 


Model 

Cp 

R2 

Adj-R2 

%RMSE  ” 

%Avg 

Err 

%RMSE 

%Max 

Err 

Stress 

Data  Set  1 (1,000  Data  Points) 

Data  Set  2 (250  Data  Points) 

Full 

Quadratic 

78.000 

0.9847 

0.9834 

-1.9226  -1.9910 

-1.6692 

-1.9128 

-4.0613 

(78  terms) 
Full 

Enhanced 

98.000 

0.9988 

0.9986 

-0.5510  -0.5794 

-0.5543 

-0.7094 

-2.2777 

(98  terms) 
Reduced 
Enhanced 

40.549 

0.9987 

0.9987 

-0.5510  -0.5622 

-0.5876 

-0.7404 

-2.1497 

(60  terms) 

example  problem,  the  reduced  response  surface  approximation  is  preferred  due  to  its  greater 
accuracy  and  smaller  number  of  parameters,  that  represents  a smoother  response  function. 

Figure  (6.11)  shows  the  correlation  between  the  finite  element  strain  values  and  the 
strain  values  obtained  from  the  response  surface  approximation.  If  there  were  100% 
correlation  between  the  finite  element  and  the  response  surface  approximation  strain  values, 
the  actual  data  points  indicated  by  the  circular  symbols  would  fall  on  the  ideal  diagonal  line. 


Figure  6.11:  Response  surface  approximation  compared  to  finite  element  strain  values 

(Data  Sets  1 and  2). 
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From  Fig.  (6.11)  it  is  clear  that  the  actual  data  points  are  close  to  the  ideal  diagonal  line, 
indicating  a high  correlation  or  accuracy  between  the  finite  element  and  the  response  surface 
approximation  strain  values. 


6.5.I.2.  Stress  distribution  response  surface  approximations  (Delamination  failurel 

The  delamination  failure  load  may  be  calculated  from  Eq.  (6.6)  and  depends  on  the 
(Tzz,  <Jzx  and  gZy  inter-laminar  stress  components  in  the  interfaces  directly  above  and  below  the 
dropped  plies.  However,  it  was  found  that  the  contribution  of  the  crZY  stress  component  was 
negligible  for  the  current  problem.  Neglecting  the  contribution  of  the  aZy  stress  component 
influenced  the  average  delamination  failure  index  value  for  all  of  the  data  points  of  both  Data 
Sets  1 and  2 by  only  0.05%,  while  the  maximum  difference  for  all  of  these  data  points  was 
only  0.6%.  Response  surface  approximations  were  thus  only  constructed  for  the  aZz  and  era- 
stress  components,  which  may  be  written  in  non-dimensional  function  form  (see  Eq.  6.10)  as: 


( h.  , . . . 

F 9 F 9 F 9 12  9 23  ? l-DroP9  l-T°P9  3_Top’  v\  J 

V -£-1  £-1  •&! 


r r r V' 

Drop"*  1 _ Top  ’ 3 _ Top  5 1 _ Bottom  ’ 3 _ Bottom  ’ 


£i_  <±n_  y V'  V’  V'  V’  - 

F 9 F 9 F 9 12  5 ^23  ? l _ Drop  5 1 _Top*  3_Top  5 l _ Bottom  * 3 _ Bottom*  j 

LL\  l j 


(6.16) 


As  was  the  case  for  the  strain  distribution  response  surface  approximation,  V3 * D depends 


linearly  on  Vt  Drop  and  is  not  required  as  a predictor  variable  (see  Section  6.5). 

For  the  stress  distributions,  the  quadratic  response  surface  approximations  were 
enhanced  by  including  all  third  order  parameters  associated  with  the  four  most  important 
predictor  variables,  which  were  identified  as  E2/Eh  GI2/Eh  V'  Drop  and  x/l.  Finally,  the 

stepwise  regression  procedure  and  the  Cp  statistic  were  used  to  reduce  the  enhanced  response 
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surface  approximations.  The  predictive  capabilities  of  all  three  response  surface 
approximations,  calculated  from  both  Data  Sets  1 and  2,  are  summarized  in  Table  (6.5). 


Table  6.5:  Stress  distributions  (a:2  and  c^)  response  surface  approximations. 


Model 

Cp 

R2 

Adj-R2 

%RMSE 

%RMSE 

PRESS 

%Avg 

Err 

%RMSE 

%Max 

Err 

°zz 

Data  Set  1 (1,200  Data  Points) 

Data  Set  2 (300  Data  Points) 

Full 

Quad 
(78  terms) 
Full 

78.000 

0.9477 

0.9441 

21.3617 

22.1134 

17.7452 

21.4668 

66.8383 

Enhanced 

98.000 

0.9893 

0.9884 

9.7205 

10.1224 

8.5347 

11.1822 

36.6367 

(98  terms) 
Reduced 

Enhanced 
(46  terms) 

14.057 

0.9892 

0.9887 

9.5937 

9.7773 

8.2640 

11.0264 

36.9316 

CTzx 

Data  Set  1 (1,200  Data  Points) 

Data  Set  2 (300  Data  Points) 

Full 

Quad 
(78  terms) 
Full 

78.000 

0.9974 

0.9973 

-2.5458 

-2.6324 

-2.4617 

-3.1395 

-7.6673 

Enhanced 

98.000 

0.9988 

0.9987 

-1.7696 

-1.8395 

-1.8273 

-2.2827 

-7.3187 

(98  terms) 
Reduced 

Enhanced 
(46  terms) 

13.942 

0.9988 

0.9987 

-1.7465 

-1.7784 

-1.8275 

-2.2718 

-6.8726 

As  for  the  yl2  response  surface  approximations,  the  enhanced  stress  distribution 
response  surface  approximations  are  more  accurate  than  the  quadratic  stress  distribution 
response  surface  approximations.  The  reduced  enhanced  response  surface  approximations  are 
slightly  more  accurate  than  the  full  enhanced  response  surface  approximations  (except  for  the 
%McixErr  value  in  the  ctzz  response  surface  approximation  and  the  %AvgErr  value  in  the  crzr 
response  surface  approximation).  Again,  the  reduced  response  surface  approximations  are 
preferred  due  to  its  greater  accuracy  and  smaller  number  of  parameters  that  represent  a 
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smoother  approximate  response  function,  and  were  used  in  the  calculation  of  the  possibility  of 
failure.  The  predictive  capabilities  of  the  reduced  response  surface  approximations  are  shown 
graphically  in  Fig.  (6.12). 


Figure  6.12:  Response  surface  approximation  compared  to  finite  element  stress  values 

(Data  Sets  1 and  2):  (a)  Gzz  distribution;  (b)  Gzx  distribution. 

From  Table  (6.5)  and  Fig.  (6.12)  it  is  clear  that  the  Gzz  stress  distribution  response 
surface  approximation  is  not  as  accurate  as  the  gZy  stress  distribution  response  surface 
approximation.  However,  the  magnitude  of  Gzz  stress  distribution  is  significantly  smaller  than 
the  Gzx  stress  distribution  and  the  value  of  the  failure  index  is  thus  dominated  by  Gzx-  The  fact 
that  the  Gzx  stress  distribution  dominates  the  delamination  failure  index  value  is  illustrated  by 
calculating  the  failure  index  from  the  two  reduced  response  surface  approximations,  with  the 
results  summarized  in  Table  (6.6)  and  shown  graphically  in  Fig.  (6.13). 


Table  6.6:  Delamination  failure  index  FDeiam  (Data  Set  2 with  300  data  points). 


%AvgErr 

%RMSE 

%MaxErr 

1.4483 

1.8173 

5.4462 
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Figure  6.13:  Delamination  failure  index  FDeiam  (calculated  from  the  reduced  crZz  and  crzr 

response  surface  approximations  for  Data  Sets  1 and  2). 


6.5.2. Possibility  of  Failure  Response  Surface  Approximations 

For  the  possibility  of  failure,  a full  cubic  polynomial  in  the  five  possibility  of  failure 
predictor  variables  (56  parameters)  was  used  as  an  initial  response  surface  approximation. 
Recall  that  the  five  possibility  of  failure  predictor  variables  are  the  weight  (fV)  and  the  level  of 
uncertainty  associated  with  the  elastic  material  properties  (EL),  the  ply  failure  strength 
allowable  ( SP ),  the  ply  lay-up  angle  (PL)  and  the  characteristic  distance  (CD).  The 
2?-Optimality  criterion  was  used  to  select  200  data  points  for  constructing  the  cubic  possibility 
of  failure  response  surface  approximations.  The  /^Optimal  data  points  were  selected  from  a 
set  of  7,776  candidate  points  that  were  created  by  considering  a six  level  factorial  design, 
consisting  of  six  equally  spaced  levels  in  each  of  the  design  variables.  Note  that  a six  level 
factorial  design  is  not  required  to  generate  candidate  points  for  a cubic  response  surface 
approximation,  but  was  used  here  since  the  small  number  of  predictor  variables  allowed  the 
use  of  six  data  points  in  each  design  variable  for  generating  a comprehensive  set  of  candidate 
points  for  selecting  the  ZTOptimal  data  points. 
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Two  hundred  data  points  are  a fairly  large  number  of  data  points  for  constructing  a 
polynomial  with  only  56  parameters.  However,  as  in  the  isotropic  example  problem  only 
data  points  with  a possibility  of  failure  not  equal  to  either  0 or  1 were  used  in  constructing  the 
possibility  of  failure  response  surface  approximations.  Since  it  is  not  possible  to  know 
beforehand  how  many  data  points  would  be  available  for  constructing  the  response  surface 
approximations,  a relatively  large  number  of  data  points  were  thus  selected.  Two  possibility 
of  failure  response  surface  approximations  were  constructed  according  to  the  two  failure 
modes.  The  resulting  predicted  possibility  of  failure  is  then  obtained  from: 


n(P_, f)(EL,SP,PL,CD,W)  = 

min( nVlf(  (EL,  SP,  PL,  CD,  W),flDelamFall  (EL,  SP,  PL,  CD,  W)) 


(6.17) 


For  the  ply  failure  possibility  of  failure,  121  of  the  original  200  data  points  had  a 
possibility  not  equal  to  either  0 or  1,  while  197  of  the  data  points  had  a delamination  failure 
possibility  of  failure  not  equal  to  either  0 or  1.  As  before,  the  general  response  surface 
approximations  were  reduced  using  the  mixed  stepwise  regression  procedure  and  the  Cp 
statistic,  with  the  predictive  capabilities  of  the  possibility  of  failure  response  surface 
approximations  summarized  in  Table  (6.7). 

As  in  the  isotropic  example  problem,  only  a single  data  set  was  used  to  construct  the 


possibility  of  failure  response  surface  approximations  and  the  %RMSEPRess  value  was  used  to 
compensate  for  the  possibly  overly  optimistic  error  estimates. 
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Table  6.7: 


Possibility  of  failure  response  surface  approximation. 


Model 

Cp  R2  Adj-R2  %RMSE 

Ply  Failure 

Third  Order  Model  (161  Data  Points) 

Full 

(56  Terms) 
Reduced 
(38  Terms) 
Delamination 
Failure 


56.000  0.9999  0.9997  0.9811  1.7165 

29.783  0.9998  0.9998  0.9313  1.3871 

Third  Order  Model  (164  Data  Points) 


Full 

(56  Terms) 
Reduced 
(40  T erms) 


56.000 

30.982 


0.9979 

0.9978 


0.9970 

0.9972 


3.1043  4.1788 

3.0139  3.8650 


CHAPTER  7 

DROPPED  PLY  COMPOSITE  LAMINATED  PLATE  OPTIMIZATION 


The  design  for  uncertainty  problem  considered  here,  involves  minimizing  a cost 
function  of  a typical  aircraft  structure,  subject  to  uncertain  problem  parameters,  for  a specified 
allowable  possibility  of  failure.  A dropped  ply  composite  laminated  plate  is  considered  as  a 
typical  structural  component  on  such  an  aircraft.  As  was  the  case  in  the  isotropic  plate 
example  problem,  it  is  assumed  that  the  structure  will  be  built  well  into  the  future  from 
materials  not  yet  available.  The  cost  function  includes  costs  associated  with  increased 
structural  weight,  costs  associated  with  material  property  characterization  and  costs  associated 
with  laying  up  the  plies  within  a specified  tolerance.  The  material  characterization  and 
manufacturing  tolerance  costs  are  applicable  to  the  entire  aircraft.  The  costs  associated  with 
increased  weight  are  thus  scaled  to  the  entire  aircraft,  by  assuming  that  the  dropped  ply 
composite  laminated  plate  is  representative  of  other  structural  components. 

Fuzzy  set  theory  is  used  to  model  the  uncertainty  using  the  limited  data  available  to 
the  designer.  Additionally,  response  surface  approximations  are  used  throughout  the  design 
process,  mainly  to  reduce  the  computational  burden  associated  with  designing  for  uncertainty, 
but  also  to  simplify  the  integration  of  the  analysis  code  with  the  optimization  algorithm. 
Response  surface  approximations  also  filter  out  numerical  noise  in  the  approximate  response 
function,  thus  allowing  the  use  of  a derivative-based  optimization  algorithm.  The  design 
optimization  trades  off  costs  associated  with  material  characterization  versus  costs  associated 
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with  increased  structural  weight.  The  optimum  results  are  presented  in  the  form  of  design 
charts,  trading  off  the  cost  against  the  allowable  possibility  of  failure. 


7.1.  General  Design  Problem 

A composite  laminated  plate  with  a change  in  thickness  across  its  width,  which  is  a 
result  of  the  internal  termination  of  plies  (see  Fig.  7.1),  is  the  present  design  problem  (more 
detail  is  presented  in  CHAPTER  6).  It  is  assumed  that  the  material  properties  and  the  applied 
load  are  uncertain,  while  manufacturing  tolerances  are  modeled  by  considering  uncertainty  in 
the  ply  lay-up  angles. 


Normal  (Z) 


Transverse  (Y) 
Longitudinal  (X) 


Continuous  Plies 
(Quasi-Isotropic) 

Drop  Plies 
(0°  Plies) 

Continuous  Plies 
(Quasi-Isotropic) 


Thick  Section  Transition  Section  Thin  Section 


Figure  7.1:  Two  step  dropped  ply  composite  laminated  plate. 


Two  failure  mechanisms  are  considered  during  the  design,  delamination  in  the 
interfaces  directly  above  and  below  the  dropped  plies,  and  ply  failure.  When  considering  both 
delamination  and  ply  failure  modes,  the  failure  load  Pf  of  the  structure  may  be  written  as 
(see  CHAPTER  6) 


Pf  -Min 


' P'y 


1.5 


(7.1) 
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where: 


p e 

p — 12 

rPly  ~ 


Yn 


P, 


Delam 


V^3  J 


9 zz 


(7.2) 


V J 


In  Eqs.  (7.1)  and  (7.2)  PP,y  denotes  the  ply  failure  load,  PDeiam  the  delamination  failure  load  and 
P the  applied  load.  Furthermore,  yn  and  en  are  the  shear  strain  and  shear  strain  allowables, 
while  &zx>  crZz , SI2  and  Z,  are  the  inter-laminar  stress  components  and  stress  allowables 
(material  properties  summarized  subsequently  in  Table  7.1),  respectively.  Failure  is  then 
defined  to  occur  when: 


P-Pf>  0 (7.3) 

The  factor  of  safety  of  1.5  applied  to  the  delamination  failure  load  in  Eq.  (7.1)  is  used  to 
account  for  the  fact  that  the  delamination  failure  load  represents  an  ultimate  load  at  which 
catastrophic  failure  occurs.  In  contrast,  the  ply  failure  load  represents  a limit  load,  indicating 
the  onset  of  failure. 

The  objective  of  the  design  problem  is  to  minimize  a cost  function  (see  Section  7.2), 
subject  to  an  allowable  possibility  of  failure.  The  results  are  presented  in  the  form  of  design 
charts,  trading  of  the  cost  against  the  allowable  possibility  of  failure.  The  problem  has  five 
design  variables,  namely  the  weight  (W)  and  the  uncertainty  associated  with  the  elastic 
material  properties  (EL),  the  ply  failure  strength  allowable  ( SP ),  the  characteristic  distance 
(CD),  and  the  ply  lay-up  angle  (PL).  The  uncertainty  in  the  delamination  strength  allowables 
(SD)  is  also  considered,  however  it  is  not  a design  variable,  but  is  assumed  to  be  twice  the 
uncertainty  associated  with  the  ply  failure  strength  allowable.  A fixed  uncertainty  associated 
with  the  applied  load  of  ± 20%  is  assumed.  The  uncertain  problem  parameters  and  the 
nominal  level  of  uncertainty  associated  with  each  parameter  are  summarized  in  Table  (7.1). 
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Table  7.1:  Uncertain  problem  parameters. 


Variable 

Nominal  Values 

Maximum 

Uncertainty 

Design  Variables 

Ei 

128,000  GPa 

±10% 

e2,e3 

11,300  GPa 

±10% 

Gi2,  G3i 

6,000  GPa 

±10% 

Elastic  Properties 

g23 

3 380  GPa 

±10% 

{EL) 

Vl2,  VI3 

0.3 

±10% 

V23 

0.35 

±10% 

z, 

52.0  MPa 

+40%t 

Strength  Properties 

s13 

92.1  MPa 

+40%* 

ei2 

0.01552 

±20% 

{SP  and  SD) 

d0 

0.25  mm 

±20% 

Characteristic  Distance 
{CD) 

p 

10  000  N 

±20% 

Variable 

Nominal  Values 

Maximum 

Uncertainty 

Design  Variables 

f) 

Default  Lay-up 

±5° 

Ply  Lay-Up 

of  sub-laminates 

m 

t The  uncertainty  associated  with  the  delamination  strength  allowables  (SDJ 
is  assumed  to  be  twice  the  uncertainty  associated  with  the  ply  strength  allowable  (SP). 


The  uncertainty  in  the  lay-up  angle  0,  influences  the  five  lamination  parameters 
K Drop*  V\_toP  > Pj  Top’  K_Boi  and  V3  Bol  in  the  reduced  stress  and  strain  distribution  response 
surface  approximations  used  to  calculated  the  failure  load  Pf.  Assuming  that  E2=E3,  Gn=Gn 
and  Vi2=V]3  (as  shown  in  Table  7.1),  the  possibility  of  failure  due  to  the  ply  and  delamination 
failure  modes  is  calculated  from  a total  of  16  uncertain  problem  parameters:  Ej,  E2,  Gi2,  G23, 
t^3,  Zr,  S;3,  e/2,  Vx  Drop,  VXTop  , V3  Top,  P\_Boi > V3_Bot  > do  and  P. 


7.2.  Cost  Function 


The  cost  function  used  in  the  present  dissertation  includes  costs  associated  with 
increased  structural  weight,  costs  associated  with  material  property  characterization  and  costs 
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associated  with  laying  up  the  plies  within  a specified  tolerance.  The  cost  function  may  be 
written  as 

Cost  = C,  Weight  + C2  Material  + C3  PlyLayUp  (7.4) 

where  the  coefficients  Ch  C2  and  C3  denotes  the  cost  coefficients  of  the  different  components 
and  will  be  discussed  in  detail  in  the  rest  of  this  section.  For  the  design  for  uncertainty 
problem,  the  design  variables  are  the  five  predictor  variables  (EL,  SP,  CD,  PL,  and  W)  of  the 
possibility  of  failure  response  surface  approximations. 

In  combining  costs  associated  with  weight,  with  costs  associated  with  characterizing 
material  properties  and  with  manufacturing  tolerances,  a problem  arises.  The  weight 
calculated  is  for  a small  structural  component,  while  the  material  property  characterization 
and  manufacturing  tolerance  costs  can  be  amortized  over  the  entire  aircraft,  in  fact,  an  entire 
fleet  of  aircraft  made  from  the  same  materials.  Here  this  problem  is  addressed  by  assuming 
that  the  dropped  ply  composite  laminated  plate  is  a representative  structural  component  and 
its  structural  weight  is  thus  scaled  to  an  entire  aircraft.  This  implicitly  assumes  that  the 
number  of  aircraft  that  will  be  built  from  the  same  material  is  approximately  equal  to  the 
number  of  material  systems  on  one  aircraft,  which  may  exaggerate  the  material 
characterization  costs. 

The  numbers  used  for  the  cost  calculation  in  the  present  dissertation  are  intended  to 
be  only  rough  estimates,  used  here  to  demonstrate  a methodology  rather  than  performing  a 
true  cost  calculation.  A realistic  cost  calculation  would  require  data  from  more  structural 
components  and  more  precise  amortization  of  the  cost  of  characterizing  material  properties  or 
reducing  manufacturing  tolerances. 
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For  the  manufacturing  and  fuel  costs  associated  with  the  structural  weight,  it  was 
assumed  that  the  failure  load  Pf  is  proportional  to  the  representative  composite  material 
weight  W of  an  entire  aircraft  as  follows: 

Pf  °c  W (7.5) 

Additionally,  if  the  representative  composite  material  weight  obtained  from  the  baseline  panel 
with  no  uncertainty  is  denoted  by  Wa  and  the  corresponding  failure  load  by  PfQ,  the  failure 
load  with  uncertainty,  Pf,  may  then  be  written  as  follows: 


(7.6) 


The  baseline  weight  Wa  of  the  structure  is  calculated  by  assuming  that  a typical  aircraft  has  a 
structural  composite  material  surface  area  of  about  1,161  m2  (12,500  ft2),  corresponding  to 
180,000  of  the  254.0  mm  by  25.46  mm  panels  considered  here.  The  weight  of  a single  panel  is 
calculated  by  multiplying  the  volume  of  the  baseline  panel  VPaneI  (VPane,=  14,651mm3)  by  the 
density  of  the  laminate  p (p=  1,589  kg/m3,  Gibson94  1994,  pp.  8 and  81).  The  baseline  weight 
Wa  is  then  obtained  as  follows: 

1T0=  180,000  (VPanel  p)= 4, 190  kg  (7.7) 

The  cost  coefficient  for  the  weight  component,  Ch  is  composed  of  the  manufacturing 
and  fuel  cost  per  kilogram  of  structural  weight,  where  the  fuel  cost  is  calculated  over  the  entire 
lifetime  of  the  aircraft.  It  was  assumed  that  a kilogram  of  structural  weight  is  worth  $550  in 
terms  of  material  and  manufacturing  cost  plus  about  0.1  kilogram  of  fuel  (or  2.94  cents)  per 
flight.  Assuming  that  the  aircraft  has  a life  span  of  25  years  and  will  complete  2,000  flights  per 
year,  one  kilogram  of  structural  weight  results  in  $1,470  of  fuel  costs  over  the  life  span  of  the 
aircraft.  The  total  manufacturing  and  fuel  cost  component  for  structural  weight  is  then: 
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C,  Weight  = 2,020  W 


(7.8) 


The  material  characterization  costs  consist  of  several  components,  corresponding  to 
the  different  material  properties.  It  was  assumed  that  the  cost  is  inversely  proportional  to  the 
level  of  uncertainty  associated  with  the  material  properties  as  follows 


C2  Material  = C2 


JO 

EL 


- + c. 


20 

SP 


- + c. 


40 

SD 


■ + cA 


20 

CD 


(7.9) 


where  the  c,  coefficients  denote  the  relative  cost  of  the  different  material  properties  and  are 
summarized  in  Table  (7.2).  Percentage  values  are  used  to  represent  the  uncertainty- 
components  EL,  SP,  SD  and  CD.  Note  that  Eq.  (7.9)  is  scaled  to  be  equal  to  C2,  when  each 
uncertainty  component  is  at  its  upper  limit  of  uncertainty. 

Table  7 .2:  Relative  cost  coefficients  for  the  material  characterization  costs. 


Property 

Value 

Cl 

V* 

c2 

V* 

c3 

1/4 

c4 

1/4 

The  overall  cost  coefficient  C2  was  determined  by  assuming  that  $2,000,000  were 
required  to  determine  the  material  properties  with  uncertainty  equal  to  the  upper  limit  shown 
in  Table  (7.1).  The  final  cost  associated  with  the  uncertainty  in  the  material  properties  may 
thus  be  written  as 


C2  Material  = 2,000,000 


10  20 
- + - 


4 EL  4 SI 


40  20  ' 

4SD  + 4CDy 


(7.10) 


The  cost  of  reducing  the  tolerance  in  the  ply  lay-up  angles  was  again  assumed  to  be 
inversely  proportional  to  the  level  of  uncertainty  associated  with  the  lay-up  angles.  For  this 
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case,  it  was  assumed  that  laying  up  the  plies  within  a tolerance  of  5°  would  correspond  to  a 
cost  of  $50,000  and  this  cost  component  may  then  be  written  as 

C3  PlyLayUp  = 50,000  (7.1 1) 

where  PL  denotes  the  uncertainty  in  the  ply  lay-up  angles  measured  in  degrees. 

7.3.  Cost  Optimization 

In  the  design  for  uncertainty  problem  considered  here,  the  cost  function  of  Section  7.2 
was  minimized,  subject  to  an  allowable  possibility  of  failure  as  a result  of  ply  and 
delamination  failure  modes.  The  possibility  of  failure  was  obtained  from  the  reduced 
possibility  of  failure  response  surface  approximations  of  CHAPTER  6.  The  optimization  was 
performed  using  the  sequential  linear  programming  subroutine  provided  with  the  commercial 
software  package,  Design  Optimization  Tools70  (DOT)  (1995).  In  order  to  avoid  local 
minima,  each  optimization  was  repeated  ten  times  from  different  starting  points  randomly 
distributed  throughout  the  design  space.  No  local  minima  were  found.  Note  that  the 
optimization  algorithm  used  here  is  different  from  that  used  in  the  isotropic  plate  example 
problem  (Microsoft  Excel  Version  8.069  1997)  and  illustrates  the  fact  that  response  surface 
approximations  allow  easy  integration  of  the  analysis  code  with  any  appropriate  optimization 
algorithm.  For  the  present  example  problem  the  Design  Optimization  Tools  (DOT)70  (1995) 
software  was  preferred  since  a large  number  of  optimizations  were  required  which  could  easily 
be  performed  from  within  a simple  FORTRAN71  (1995)  program. 

7.3.1. Optimization  Problem 

The  optimization  problem  for  minimizing  the  cost  function  subject  to  a allowable 
possibility  of  failure,  may  be  defined  as: 
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Minimize: 


Cost  = 2,020  W + 2,000,000 


10  20  40  20  "l  fnA„J  5 ^ 

■ H 1 1 I + 50,000 

PL 


4 EL  4 SP  4 SD  4 CD  ) 


(7.12) 


Subject  to: 


nfrr/)(EL,SP,PL,CIXW) 


n 


-1<0 


allowable 


(7.13) 


Where  n(P_P/)  denotes  the  approximate  possibility  of  failure  as  obtained  from  the  reduced 
response  surface  approximation  of  CHAPTER  6 and  n „ u denotes  the  allowable 
possibility  of  failure.  The  value  of  Llallmvable  was  varied  from  0 to  1 to  study  the  influence  of 
the  allowable  possibility  of  failure  on  the  cost  function  and  its  individual  cost  components. 

7.3.2. Results 

First,  the  optimization  was  performed  for  a single  n iff  M ifJj  value  equal  to  0.2.  This 
problem  was  solved  to  evaluate  the  methodology  as  well  as  the  accuracy  of  the 
approximations,  and  the  detailed  results  are  summarized  in  Table  (7.3). 

For  an  allowable  possibility  of  failure  value  equal  to  0.2,  the  optimum  cost  function  is 
equal  to  $14,218  M.  The  total  cost  is  dominated  by  the  structural  weight  cost  component 
($10,435  M or  73.4%).  The  structural  weight  is  increased  from  the  baseline  value  of  4,190  kg 
to  5,166  kg  to  increase  the  failure  load  of  the  structure,  thus  reducing  the  possibility  of  failure. 
Note  that  the  optimizer  limited  the  increase  in  the  structural  weight,  by  decreasing  the 
uncertainty  in  the  problem  parameters,  particularly  the  ply  and  delamination  failure  strength 
allowables.  This  decrease  in  the  uncertainty  associated  with  the  problem  parameters  increased 
the  material  characterization  and  manufacturing  tolerance  costs  from  $2.05  M to  $3.78  M.  If 
the  material  characterization  and  design  tolerance  costs  were  not  allowed  to  increase  above 
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Table  7.3:  Fuzzy  set  based  optimum  results  for  UallowMe  =0.2. 


Variable Value  Cost 


EL 

7.18  % 

$0.6965  M 

SP 

9.73  % 

$1.0280  M 

SD 

19.5  % 

$1.0280  M 

CD 

14.5  % 

$0.6906  M 

PL 

0.73° 

$0.3402  M 

Load 

20.0  % 

- 

Weight 

5,165.7  kg 

$10,435  M 

Cost 

-- 

$14,218  M 

nPiv 

0.0000 

(Inactive) 

(0.0000) 

^Delam 

0.2001 

(Active) 

(0.1856) 

Values  in  parenthesis  are  predicted  bp  possibility  of  failure  response  surface  approximation. 

$2.05  M,  the  structural  weight  would  have  increased  to  6,915  kg  and  the  overall  cost  to 
$16.02  M.  The  possibility  of  failure  value  calculated  from  the  vertex  method  and  the  value 
predicted  from  the  reduced  possibility  of  failure  response  surface  approximations  (values  in 
parenthesis)  are  both  shown  in  Table  (7.3),  with  the  maximum  difference  equal  to  7.2  %. 

Next,  a design  chart  of  the  cost  compared  to  the  allowable  possibility  of  failure  was 
generated.  The  allowable  possibility  of  failure  was  varied  from  0 to  1 with  a constant  interval 
of  0.05  and  the  results  are  shown  graphically  in  Fig  (7.2).  From  Fig.  (7.2)  it  is  clear  that  a 
decrease  in  the  allowable  possibility  of  failure  results  in  an  increase  in  the  cost  function.  The 
relationship  between  the  allowable  possibility  of  failure  and  the  cost  function  is  almost  linear 
with  a change  in  the  slope  at  an  allowable  possibility  of  failure  equal  to  roughly  0.8.  The 
practical  application  region  of  Fig.  (7.2)  is  for  an  allowable  possibility  of  failure  between  0 and 
0.4.  Increasing  the  allowable  possibility  of  failure  above  0.4  would  be  impractical  due  to  a 
high  possibility  of  failure. 
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Figure  7.2:  Cost  compared  to  the  allowable  possibility  of  failure. 

The  change  in  the  slope  of  the  cost  versus  the  allowable  possibility  of  failure  design 
chart  at  an  allowable  possibility  of  failure  roughly  equal  to  0.8  corresponds  to  a change  in  the 
active  failure  mode,  as  shown  in  Fig.  (7.3).  The  normalized  constraint  functions  of  Fig.  (7.3) 
are  inactive  when  negative,  active  when  equal  to  zero  and  violated  when  positive. 
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Figure  7.3:  Normalized  constraint  functions  compared  to  the  allowable  possibility  of 

failure  (active  when  normalized  constraint  function  is  equal  to  0). 

Finally,  the  cost  associated  with  each  of  the  design  variables  is  shown  as  a function  of 
the  allowable  possibility  of  failure  in  Fig.  (7.4).  As  expected,  the  cost  associated  with  the 
structural  weight  dominates  the  overall  cost  component.  The  cost  associated  with  the 
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Figure  7.4:  Individual  cost  components  compared  to  the  allowable  possibility  of  failure. 

strength  allowables  dominates  the  material  characterization  and  manufacturing  tolerance 
costs.  Of  all  the  uncertainties  associated  with  the  uncertain  problem  parameters,  the 
uncertainty  associated  with  the  strength  allowables  changes  the  most  with  a change  in  the 
allowable  possibility  of  lailure.  Reducing  the  uncertainty  in  the  elastic  material  properties  and 
the  characteristic  distance  costs  about  the  same,  while  reducing  the  uncertainty  in  the  ply  lay- 
up angles  is  fairly  inexpensive.  From  Fig.  (7.4)  it  is  clear  that  the  optimizer  increases  the 
structural  weight  with  a decrease  in  the  allowable  possibility  of  failure  in  order  to  increase  the 
failure  load  of  the  structure. 

As  noted  previously,  the  material  characterization  costs  can  be  amortized  over  many 
aircraft  that  are  built  from  the  same  materials.  In  the  present  dissertation,  it  was  implicitly 
assumed  that  the  number  of  aircraft  that  will  be  built  from  the  same  materials  is  equal  to  the 
number  of  material  systems  for  a single  aircraft.  Since  this  assumption  may  exaggerate  the 
material  characterization  costs,  a second  set  of  optimizations  was  performed  to  demonstrate 
the  cost  calculation  for  the  case  where  more  aircraft  will  be  built  from  the  same  materials.  In 
this  case  it  was  assumed  that  the  number  of  aircraft  that  will  be  built  from  the  same  materials 
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is  equal  to  ten  times  the  number  of  material  systems  for  a single  aircraft,  thus  reducing  the 
material  characterization  costs  by  a factor  of  ten.  The  cost  coefficient  C2  of  Eq.  (7.4)  is  thus 
reduced  from  2,000,000  to  200,000  and  the  material  characterization  costs  are  thus  defined  as: 


C2  Material  = 200,000 


( 10  20  40  20  N 


- + - 


- + - 


(7.14) 


4 EL  ASP  4 SD  4 CD 

Detailed  results  for  the  cost  function  with  reduced  material  characterization  costs 
(see  Eq.  7.14)  are  summarized  in  Table  (7.4)  for  the  case  where  the  allowable  possibility  of 
failure  is  equal  to  the  0.2.  As  expected,  the  reduced  material  characterization  costs  resulted  in 
a reduction  in  the  level  of  uncertainty  associated  with  the  material  properties.  This  reduction 
in  the  level  of  uncertainty  associated  with  the  material  properties,  resulted  in  a decrease  in 
both  the  structural  weight  cost  component  and  the  manufacturing  tolerance  cost  component. 
The  costs  associated  with  the  structural  weight  were  decreased  from  $10,435  M to  $8,798  M, 
the  costs  associated  with  material  characterization  were  reduced  from  $3,443  M to  $0,879  M, 


Table  7.4:  Fuzzy  set  based  optimum  results  for  11,  „ =0.2. 


Variable 

Value 

Cost 

EL 

2.80  % 

$0.1783  M 

SP 

3.85  % 

$0.2597  M 

SD 

7.70  % 

$0.2597  M 

CD 

5.51  % 

$0.1815  M 

PL 

0.91° 

$0.2760  M 

Load 

20.0  % 

- 

Weight 

4,355.3  kg 

$8.7976  M 

Cost 

- 

$9.9530  M 

nply 

0.2000 

(Active) 

(0.1962) 

“Ddam 

0.2000 

(Active) 

(0.1818) 

Values  in  parenthesis  are  predicted  by  possibility  of  failure  response  surface  approximation. 
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and  the  costs  associated  with  the  manufacturing  tolerances  were  reduced  from  $0,340  M to 
$0,276  M.  The  overall  cost  was  reduced  from  $14,218  M to  $9,953  M.  The  greatest  reduction 
in  uncertainty  was  observed  for  the  uncertainty  associated  with  the  strength  allowables  and 
the  characteristic  distance.  The  reduction  in  the  uncertainty  associated  with  the  material 
properties  resulted  in  a slight  increase  in  the  tolerance  of  the  ply  lay-up  angle  and  a reduction 
in  the  structural  weight.  The  tolerance  in  the  ply  lay-up  angle  was  increased  from  0.73°  to 
0.91°  and  the  structural  weight  was  reduced  from  5,166  kg  to  4,355  kg.  The  possibility  of 
failure  value  calculated  from  the  vertex  method  and  the  value  predicted  from  the  reduced 
possibility  of  failure  response  surface  approximations  (values  in  parenthesis)  are  both  shown  in 
Table  (7.4).  For  ply  failure,  the  values  of  the  possibility  of  failure  obtained  from  the  vertex 
method  and  the  reduced  response  surface  approximation  respectively  differs  by  1.9%,  while 
the  difference  is  equal  to  9.1%  for  delamination  failure. 

The  cost  function  is  shown  graphically  as  a function  of  the  allowable  possibility  of 
failure  in  Fig.  (7.5).  The  allowable  possibility  of  failure  was  varied  from  0 to  1,  with  a 
constant  interval  of  0.05.  As  before,  a decrease  in  the  allowable  possibility  of  failure  results  in 


Allowable  Possibility  of  Failure 


Figure  7.5:  Cost  compared  to  the  allowable  possibility  of  failure. 
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an  increase  in  the  cost  function,  with  an  almost  linear  relationship  between  the  allowable 
possibility  of  failure  and  the  cost  function.  The  maximum  value  of  the  cost  function  is, 
however,  reduced  from  roughly  $15.3  M to  $10.5  M. 

The  normalized  constraint  functions  are  shown  graphically  in  Fig.  (7.6)  as  a function 
of  the  allowable  possibility  of  failure.  In  this  case,  both  constraints  are  active  over  a wide 
range  of  allowable  possibility  of  failure  values.  However,  for  small  values  of  the  allowable 
possibility  of  failure,  only  the  delamination  failure  constraint  is  active,  while  only  the  ply 
failure  constraint  is  active  for  large  values  of  the  allowable  possibility  of  failure. 


Allowable  Possibility  of  Failure 

Figure  7.6:  Normalized  constraint  functions  compared  to  the  allowable  possibility  of 

failure  (active  when  normalized  constraint  function  is  equal  to  0). 

Finally,  the  individual  cost  components  of  the  objective  function  (the  cost  function) 
are  shown  graphically  in  Fig.  (7.7)  as  a function  of  the  allowable  possibility  of  failure.  The 
individual  cost  components  are  represented  by  curves  that  are  similar  to  those  of  the  initial 
example  problem.  As  expected,  the  absolute  values  of  the  material  characterization  costs  are 
much  lower,  while  the  structural  weight  cost  component  is  even  more  dominant  in  terms  of 
the  overall  cost.  Note,  however,  that  the  absolute  value  of  the  structural  weight  cost 
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Allowable  Possibility  of  Failure 


Figure  7.7:  Individual  cost  components  compared  to  the  allowable  possibility  of  failure. 

component  is  lower  than  in  the  initial  problem,  due  to  more  accurate  material  property 
characterization. 

7.4.  Reduction  in  Computational  Cost 

The  reduction  in  computational  cost  is  mainly  associated  with  the  response  surface 
approximations  that  were  constructed  for  the  stress  and  strain  distributions.  These  response 
surface  approximations  of  CHAPTER  6 were  used  to  replace  computationally  expensive  finite 
element  analyses  during  the  evaluation  of  the  possibility  of  failure.  Compared  to  the  isotropic 
plate  example  problem,  this  process  resulted  in  even  more  substantial  reductions  in 
computational  cost. 

The  reductions  in  computational  cost  are  illustrated  for  a single  two  level  optimization 
problem.  The  evaluation  of  the  possibility  of  failure  for  a single  a level  cut  requires 
210*20*20+211*20*20=  1,228,800  finite  element  analyses  when  no  response  surface 
approximations  are  used.  Recall  that  there  are  a total  of  16  uncertain  problem  parameters, 
14  associated  with  the  ply  failure  mode  and  15  associated  with  the  delamination  failure  mode. 
Additionally,  it  was  necessary  to  consider  all  boundary  points  (20  data  points)  for  each  of  the 
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K _Top  and  K Top  ^ weh  as  the  V1Bot  and  Vj_Bol  design  spaces,  respectively  (see  CHAPTER  6). 

When  not  using  response  surface  approximations,  an  estimate  of  the  required  number  of  finite 
element  analyses  to  perform  a single  optimization  is  obtained  from  the  product  of  the 
following  four  numbers: 

Average  number  of  design  optimization  iterations:  10 

Average  number  of  n(P_Py)  evaluations  per  iteration:  10 

Average  number  of  a level  cut  evaluations  per  evaluation:  5 

Number  of  finite  element  analyses  per  a level  cut  evaluation  of  11^.^:  1,228,800 
Total  number  of  finite  element  analyses  required  per  optimization:  614,400,000 

In  contrast,  the  stress  distribution  and  buckling  load  response  surface  approximations 
were  constructed,  and  their  predictive  capabilities  evaluated,  from  a total  of  only  250  finite 
element  analyses.  These  observations  indicates  that  response  surface  approximations  is  an 
effective  approach  for  reducing  the  computational  cost  associated  with  performing  fuzzy  set 
based  designing  for  uncertainty. 

As  in  the  isotropic  plate  example  problem,  the  use  of  response  surface  approximations 
allowed  multiple  optimizations  to  be  conducted  at  little  additional  computational  cost.  It  was 
thus  possible  to  construct  and  investigate  design  charts  of  the  optimum  results,  based  on 
multiple  optimizations. 


CHAPTER  8 

CONCLUDING  REMARKS 


The  use  of  response  surface  approximations  in  expensive  structural  optimization 
problems  was  investigated  in  this  dissertation.  Designing  for  uncertainty  of  a stepped  plate 
was  considered  as  a representative  example  of  an  expensive  structural  optimization  problem. 
Both  an  isotropic  and  a composite  laminated  structure  were  considered.  For  the  composite 
laminated  structure,  the  change  in  thickness  is  a result  of  the  internal  termination  of  plies.  It 
was  assumed  that  the  structures  would  be  built  well  into  the  future,  from  materials  not  yet 
available.  Little  information  is  thus  available  on  the  uncertainty,  and  thus  fuzzy  set  theory 
was  used  to  model  the  uncertainty  in  the  material  properties,  the  geometry  and  the  applied 
load.  Designing  for  uncertainty,  where  the  uncertainty  is  modeled  by  fuzzy  set  theory,  is  a 
good  example  of  an  expensive  structural  optimization  problem,  since  it  is  a two  level 
optimization  problem. 

When  using  an  approximate  response  function  in  optimization,  the  accuracy  of  the 
response  function  is  of  the  utmost  importance,  since  optimization  procedures  have  the  general 
tendency  of  exploiting  weaknesses  in  the  response  function  formulation.  In  the  present 
dissertation,  highly  accurate  response  surface  approximations  were  obtained  by  using  a variety 
of  mathematical  and  statistical  tools.  Dimensional  analysis  was  used  to  identify  the  intrinsic 
variables  of  the  response  functions.  Using  these  variables  to  construct  the  response  surface 
approximations  improved  the  accuracy  while  reducing  the  number  of  predictor  variables  of 
the  resulting  response  surface  approximations.  Additionally,  higher  order  (cubic  and  quartic 
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instead  of  the  more  traditional  quadratic)  polynomials  were  used  as  initial  response  models, 
which  resulted  in  more  accurate  response  surface  approximations  that  were  valid  over  the 
entire  design  space.  Mixed  stepwise  regression  and  the  Cp  statistic  were  used  to  eliminate 
insignificant  parameters  from  these  initial  models,  greatly  reducing  the  number  of  parameters 
in  the  final  response  surface  approximations.  For  the  composite  laminated  plate,  where  it  was 
not  feasible  to  evaluate  a large  number  of  data  points  for  constructing  the  response  surface 
approximations,  statistical  design  of  experiments  in  the  form  of  the  ZXOptimality  criterion 
was  used  to  find  a small  set  of  data  points  while  maintaining  accuracy.  Finally,  a detailed  error 
analysis  was  performed  to  ensure  that  the  response  surface  approximations  were  indeed  of 
high  quality. 

Two  sets  of  response  surface  approximations  were  constructed  during  the  design 
process.  The  first  set  was  used  to  calculate  the  possibility  of  failure,  which  is  an  expensive 
global  optimization  problem.  In  this  case,  the  response  surface  approximations  were  mainly 
used  to  reduce  the  computational  burden  associated  with  calculating  the  possibility  of  failure. 
The  second  set  of  response  surface  approximations  was  constructed  to  represent  the  possibility 
of  failure  in  terms  of  the  overall  design  variables  of  the  problem  and  was  mainly  used  to 
simplify  the  integration  of  the  analysis  code  and  the  optimization  algorithm.  Overall,  the 
response  surface  approximations  resulted  in  substantial  reductions  in  computational  cost  since: 

1.  It  was  cheaper  to  construct  and  to  evaluate  the  predictive  capabilities  of  the 
approximations  than  to  perform  a single  two  level  optimization  problem. 

2.  The  computational  cost  was  shifted  from  the  optimization  problem  to  that  of 
constructing  the  approximations,  thus  allowing  multiple  optimizations  to  be 
conducted  at  little  additional  computational  cost. 
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The  reduction  in  computation  cost,  especially  shifting  the  cost  from  the  optimization 
problem  to  that  of  constructing  the  approximations,  allowed  the  construction  of  weight 
versus  cost  design  charts  based  on  the  results  obtained  from  multiple  optimizations.  These 
design  charts  provide  a valuable  tool  during  the  initial  design  stages.  In  the  present 
dissertation,  these  weight  versus  cost  design  charts  would  have  been  very  costly  to  construct 
without  the  use  of  response  surface  approximations. 

The  main  problem  associated  with  the  use  of  response  surface  approximations  in 
practical  problems  is  that  the  cost  of  constructing  the  approximations  is  subject  to  the  curse  of 
dimensionality.  The  curse  of  dimensionality  implies  that  the  number  of  data  points  required 
to  construct  a response  surface  approximation  is  greatly  increased  with  an  increase  in  the 
order  of  the  response  surface  approximation  and/or  the  number  of  design  variables.  Response 
surface  approximations  are  thus  not  applicable  to  all  problems,  and  the  practical  application  of 
response  surface  approximations  is  limited  to  problems  with  less  than  25  to  30  design 
variables. 

In  general,  it  has  been  illustrated  that  response  surface  approximations  provide  an 
effective  approach  for  solving  expensive  structural  optimization  problems  with  a relatively 
small  number  of  design  variables.  Response  surface  approximations  resulted  in  substantial 
reductions  in  computational  cost,  greatly  simplified  the  integration  of  the  analysis  codes  and 
optimization  algorithms,  and  allowed  the  construction  of  weight  versus  cost  design  charts. 


APPENDIX  A:  DROPPED  PLY  COMPOSITE  LAMINATED  PLATE  PROBLEM 

EVOLUTION 


A.l.  Introduction 

The  internal  termination  of  plies  (referred  to  as  dropped  plies)  is  an  effective  method 
for  reducing  the  thickness  of  composite  laminated  structures.  Unfortunately,  the  dropped 
plies  result  in  severe  stress  concentrations  due  to  the  stiffness  and  thickness  discontinuities 
associated  with  the  internal  termination  of  the  plies.  These  stress  concentrations  may  greatly 
reduce  the  strength  of  the  structure,  thus  offsetting  the  advantages  of  reducing  the  thickness 
by  terminating  plies.  Various  analytical,  numerical  and  experimental  studies  (Refs.  79  to  90) 
indicated  that  delamination  is  a common  mode  of  failure  for  dropped  ply  composite  laminated 
structures. 

The  final  transition  zone  configuration  of  the  dropped  ply  composite  laminated  plate 
considered  in  this  dissertation,  was  obtained  by  performing  a numerical  investigation  using 
detailed  finite  element  analyses.  Additionally,  manufacturing  constraints  were  also 
considered.  The  goal  of  the  numerical  investigation  was  to  study  the  influence  of  the 
transition  zone  configuration  on  the  strength  of  the  dropped  ply  composite  laminated  plate. 
During  the  numerical  investigation,  a problem  similar  to  that  studied  by  Curry,  Johnson  and 
Starnes82  (1992)  and  Davila  and  Johnson87  (1993)  was  used  as  a baseline  problem.  This  problem 
is  discussed  in  detail  in  CHAPTER  6 and  is  shown  schematically  in  Fig.  A.l.  For  the  baseline 
problem,  four  0°  plies  were  dropped  with  quasi-iso  tropic  sub-laminates  above  and  below  the 
dropped  plies. 
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Figure  A.l:  Baseline  dropped  ply  composite  laminated  plate. 


A.2.  Numerical  Investigation 

Based  on  the  experimental  and  numerical  work  of  Curry,  Johnson  and  Starnes82  (1992) 
and  the  numerical  work  of  Davila  and  Johnson87  (1993),  the  numerical  investigation  presented 
here  considered  delamination  failure  that  occurs  in  the  interfaces  directly  above  and  below  the 
dropped  plies  as  the  critical  failure  mode.  A quadratic  failure  index  was  used  to  predict 
delamination  failure  initiation  (see  CHAPTER  6 for  more  detail),  which  may  be  summarized 
as  follows. 
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This  failure  index  is  evaluated  at  a characteristic  distance  away  from  the  ply  drop.  Based  on 
the  results  of  Davila  and  Johnson87  (1993),  a characteristic  distance  of  0.25  mm  was  assumed  as 
a characteristic  distance  throughout  this  dissertation. 

The  influence  of  the  drop  zone  length  (shown  in  Fig.  A.2)  on  the  delamination 
strength  of  the  dropped  ply  composite  laminated  plate  was  investigated  first.  The  drop  zone 
length  was  varied  between  0.0887  mm  and  2.850  mm,  changing  the  drop  zone  angle  (shown  in 
Fig.  A.2)  from  10°  to  80°. 
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Figure  A.2:  Transition  zone  parameters  for  the  baseline  dropped  ply  composite  laminated 

plate. 

The  delamination  failure  index  compared  to  the  drop  zone  angle  is  shown  graphically 
in  Fig.  (A.3),  where  a smaller  drop  zone  angle  corresponds  to  a shorter  drop  zone  length. 


Figure  A.3:  Delamination  failure  index  compared  to  the  drop  zone  angle  for  the  baseline 

dropped  ply  composite  laminated  plate. 

Figure  (A.3)  indicates  that  an  increase  in  the  drop  zone  length,  corresponding  to  an 
increase  in  the  drop  zone  angle,  results  in  a decrease  in  the  delamination  strength  of  the 
composite  laminated  plate.  The  delamination  failure  index  is  thus  dominated  by  the  severity 
of  the  stiffness  discontinuity  that  occurs  at  the  ply  drop,  rather  than  the  severity  of  the 
thickness  discontinuity  which  becomes  less  severe  for  a longer  drop  zone  length. 

The  results  of  Fig.  (A.3)  were  validated  by  studying  the  accuracy  of  the  finite  element 
model  and  the  validity  of  using  equivalent  orthotropic  material  properties  for  each  sub- 
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laminate.  The  accuracy  of  the  finite  element  model  was  investigated  by  refining  the  finite 
element  mesh  in  order  to  evaluate  the  convergence  of  the  model.  The  validity  of  the 
equivalent  orthotropic  material  properties  assumption  was  investigated  by  modeling  each  ply 
of  each  sub-laminate  as  an  individual  orthotropic  layer,  orientated  at  a prescribed  angle  with 
respect  to  the  x-axis.  The  results  obtained  from  these  investigations  indicated  that:  (1)  the 
original  finite  element  model  was  not  converged;  and  (2)  the  use  of  equivalent  orthotropic 
material  properties  have  an  insignificant  influence  on  the  delamination  failure  index.  Based  on 
these  results,  the  finite  element  model  was  refined  and  the  refined  model  was  used  throughout 
the  rest  of  this  dissertation  (see  CHAPTER  6 for  more  detail  about  the  refined  model). 
However,  after  the  finite  element  model  was  refined,  the  delamination  strength  of  the 
dropped  ply  composite  laminated  plate  still  increased  with  a decrease  in  the  drop  zone  length. 

Based  on  the  validation  of  the  results  shown  in  Fig.  (A.3),  the  following  hypothesis  has 
been  formulated  to  explain  the  results: 

The  delamination  failure  index  value  of  a dropped  pip  composite  laminated 
plate  depends  on  the  severity  of  the  change  in  stiffness  between  the  thick  and  the  thin 
sections  ol  the  composite  laminated  plate. 

This  hypothesis  implies  that,  when  a resin  pocket  is  present,  a shorter  drop  zone 
length  results  in  a lower  delamination  failure  index,  since  the  severity  of  the  stiffness 
discontinuity  is  decreased.  The  resin  pocket  has  a much  lower  stiffness  than  the  surrounding 
sub-laminates  and  a longer  drop  zone  length  increases  the  compliance  of  the  resin  pocket.  The 
severity  of  the  stiffness  discontinuity  is  thus  increased  with  an  increase  in  the  drop  zone 
length. 
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The  hypothesis  was  evaluated  by  considering  a number  of  alternative  transition  zone 
configurations,  thus  changing  the  stiffness  of  the  transition  zone  and  the  severity  of  the 
stiffness  discontinuity.  The  stiffness  of  the  transition  zone  was  increased  by  reducing  the  size 
of  the  resin  pocket  as  well  as  by  reinforcing  the  resin  pocket  with  90°  fibers,  measured  relative 
to  the  x-axis.  The  size  of  the  resin  pocket  was  reduced  by  dropping  the  0°  plies  in  two  steps, 
which  resulted  in  two  smaller  resin  pockets,  and  by  tapering  the  dropped  plies  in  the  drop 
zone,  thus  completely  eliminating  the  resin  pocket  (see  Fig.  A.4). 


(a)  (b) 

Figure  A.4:  Modified  transition  zone  configurations:  (a)  two  step  configuration; 

(b)  tapered  dropped  ply  configuration. 


The  delamination  failure  index  for  both  a short  drop  zone  length  (i.e.,  0.0887  mm  or 
10°  drop  zone  angle)  and  a long  drop  zone  length  (i.e.,  1.273  mm  or  68.4°  drop  zone  angle)  are 
shown  graphically  in  Fig.  (A. 5)  for  the  baseline,  the  two  step,  and  the  tapered  configurations. 
For  the  baseline  and  the  two  step  configurations,  two  cases  were  considered.  In  the  first  case, 
the  resin  pocket  consisted  of  neat  resin,  and  in  the  second  case  the  resin  pocket  was  reinforced 
with  90°  fibers. 

From  Fig.  (A.5)  it  is  clear  that  as  the  stiffness  of  the  transition  zone  is  increased,  the 
delamination  failure  index  decreases.  When  the  resin  pocket  is  completely  eliminated  (tapered 
configuration)  the  results  of  Fig.  (A.3)  are  reversed,  since  a longer  transition  zone  increases  the 
delamination  strength  of  the  structure.  Additionally,  the  tapered  configuration  resulted  in  the 
lowest  value  of  the  delamination  failure  index,  corresponding  to  a significant  increase  in  the 
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Figure  A.5:  Delamination  failure  index  for  different  transition  zone  configurations. 

delamination  strength  of  the  structure.  Unfortunately,  high  manufacturing  costs  limit  the 
usefulness  of  the  tapered  configuration.  The  two  step  configuration  with  90°  fibers 
reinforcing  the  resin  pockets  is  the  best  practical  configuration.  For  all  cases,  except  the 
tapered  configuration,  it  was  found  that  a short  transition  zone  resulted  in  the  lowest  value  of 
the  delamination  failure  index.  However,  a shorter  transition  zone  increases  the  severity  of 
the  thickness  discontinuity.  Since  a shorter  transition  zone  results  in  lower  values  of  the 
delamination  failure  index,  even  though  the  severity  of  the  thickness  discontinuity  is 
increased,  it  may  be  concluded  that  the  thickness  discontinuity  has  a small  influence  on  the 
delamination  failure  index  value. 

Based  on  the  fact  that  a short  transition  zone  was  the  practical  configuration  with  the 
lowest  delamination  failure  index,  it  was  decided  to  investigate  ply  failure  as  a second  failure 
mechanism.  Ply  failure  occurs  due  to  high  in-plane  stresses  resulting  from  the  thickness 
discontinuity,  but  is  also  influenced  by  the  internal  stiffness  discontinuity.  Using  the 
maximum  strain  failure  criterion,  the  topmost  ply  of  the  thinnest  section  of  the  composite 
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laminated  plate  was  identified  as  the  most  critical  ply  (see  CHAPTER  6 for  more  detail).  In 
all  cases  the  in-plane  shear  strain  was  the  most  critical  strain  component,  with  a singularity 
occurring  at  the  re-entrant  corner  of  the  composite  laminated  plate.  Based  on  the  maximum 
strain  failure  criterion,  the  ply  failure  index  is  defined  as: 


F = 

1 Ply 


y 12 


(A.2) 


As  was  the  case  for  the  delamination  failure  index  of  Eq.  (A.l),  the  ply  failure  index  is 
also  calculated  at  a characteristic  distance  away  from  the  singularity,  in  this  case  the  re-entrant 
corner.  Throughout  this  dissertation,  the  same  characteristic  distance  (i.e.,  0.25  mm)  was 
assumed  for  both  the  delamination  and  ply  failure  modes.  The  ply  failure  index  value  is 
shown  graphically  in  Fig.  (A. 6),  for  the  five  transition  zone  configurations  shown  in 
Fig.  (A.5). 

From  Fig.  (A.6)  it  is  observed  that  increasing  the  stiffness  of  the  transition  zone  has  a 
much  smaller  influence  on  the  ply  failure  index  value  as  compared  to  the  influence  on  the 


Figure  A.6:  Ply  failure  index  for  different  transition  zone  configurations. 
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delamination  failure  index  value  (see  Fig.  A.5).  The  increased  stiffness  of  the  transition  zone 
has  almost  no  influence  on  the  ply  failure  index  value  for  the  short  drop  zone  length  cases,  but 
does  slightly  improve  the  ply  failure  index  value  for  the  long  drop  zone  length  cases.  As  was 
the  case  for  the  delamination  failure  index,  a short  drop  zone  length  is  the  practical 
configuration  with  the  lowest  ply  failure  index  value  (see  Fig.  A.6).  For  all  the  practical 
configurations,  the  severity  of  the  stiffness  discontinuity  dominates  the  ply  failure  index  value, 
since  the  severity  of  the  thickness  discontinuity  increases  with  a decrease  in  the  drop  zone 
length,  while  the  severity  of  the  stiffness  discontinuity  decreases  with  a decrease  in  the  drop 
zone  length. 

Figure  (A. 7)  shows  the  critical  failure  indices  when  considering  both  delamination  and 
ply  failure,  for  the  five  cases  shown  in  Figs.  (A.5)  and  (A.6).  In  order  to  make  a direct 
comparison  between  the  delamination  and  the  ply  failure  indices,  it  is  important  to  remember 
that  the  delamination  failure  load  represents  an  ultimate  load  condition,  at  which  catastrophic 
failure  occurs.  In  contrast,  the  ply  failure  load  represents  a limit  load  condition  that  indicates 
the  onset  of  failure.  Within  the  aircraft  industry  the  ultimate  load  is  typically  assumed  to  be 
1.5  times  larger  than  the  limit  load.  This  ratio  of  1.5  was  used  as  a factor  of  safety  in 
evaluating  the  delamination  failure  criterion,  while  a factor  of  safety  of  1.0  was  used  in 
evaluating  the  ply  failure  criterion.  Note  that  the  factor  of  safety  of  1.5  was  incorporated  in 
the  results  of  Figs.  (A.5)  and  (A. 7),  but  not  in  the  initial  results  of  Fig.  (A.3). 

In  Fig.  (A. 7)  the  ply  failure  mode  is  critical  for  all  short  drop  zone  length  cases  (except 
for  the  baseline  configuration),  while  the  delamination  failure  mode  is  critical  for  all  long  drop 
zone  length  cases  (except  for  the  tapered  configuration).  Additionally,  the  two  step 
configuration,  where  the  resin  pockets  are  reinforced  with  90°  fibers,  emerged  as  the  best 
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Figure  A .7:  Critical  failure  index  for  different  transition  zone  configurations. 

practical  configuration.  Figure  (A. 7)  indicates  that  for  the  two  step  configuration  with 
reinforced  resin  pockets,  a short  transition  zone  results  in  the  lowest  critical  failure  index 
value.  However,  the  difference  in  the  failure  indices  of  the  short  and  long  drop  zone  lengths  is 
small  and,  due  to  manufacturing  constraints,  it  was  decided  to  use  a drop  zone  length  equal  to 
eight  ply  thicknesses  (i.e.,  1.0056  mm  or  63.4°  drop  zone  angle). 

This  configuration  was  further  investigated,  by  considering  the  influence  of  separating 
the  two  ply  drops  on  the  strength  of  the  composite  laminated  plate.  Figure  (A.8)  shows  the 
separation  in  the  two  ply  drops  schematically.  For  this  part  of  the  investigation,  the 
separation  distance  was  varied  between  0 mm  and  4.5  mm  (i.e.,  between  0 ply  thicknesses  and 
35  ply  thicknesses)  and  the  results  are  shown  graphically  in  Fig.  (A.9). 


Separation  Distance 


Figure  A.8:  Two  step  configuration  with  separated  ply  drops. 
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Figure  A.9:  Failure  index  compared  to  the  separation  distance  for  the  two  step 

configuration  with  reinforced  resin  pockets. 

Figure  (A.9)  indicates  that  the  separation  distance  has  a much  higher  influence  on  the 
delamination  failure  index  value  than  on  the  ply  failure  index  value.  For  the  delamination 
failure  index,  there  is  a 32.8%  change  in  the  failure  index  values  when  changing  the  separation 
distance  from  0 to  35  ply  thicknesses,  while  there  is  only  a 5.8%  change  in  the  values  of  the 
ply  failure  index.  Increasing  the  separation  distance  above  roughly  15  ply  thicknesses  has  very 
little  influence  on  the  failure  indices,  since  the  two  stress  singularities  corresponding  to  the 
two  ply  drops  become  separated  and  do  not  influence  each  other.  However,  due  to 
manufacturing  constraints,  a separation  distance  of  20  ply  thicknesses  was  used  in  the  final 
configuration  considered  as  an  example  problem  in  this  dissertation.  This  final  two  step 
dropped  ply  composite  laminated  plate  configuration,  where  the  ply  drops  are  separated  by 
20  ply  thicknesses  and  the  resin  pockets  are  reinforced  with  90°  fibers,  is  shown  schematically 
in  Fig.  (A.  10). 


Ply  Failure 

i 1 1 1 1 i 1 1 1 1 i 1 1 v- 1 i 1 1 1 1 i 


137 


Normal  (Z) 


Continuous  Plies 
(Quasi-Isotropic) 

Drop  Plies 
(0°  Plies) 

Continuous  Plies 
(Quasi-Isotropic) 


Thick  Section 


Transverse  (Y) 

Longitudinal  (X) 

Transition  Section 


Thin  Section 


Reinforced  Resin  Pockets 
First  Ply  Drop  Second  Ply  Drop 


Figure  A.10:  Two  step  configuration  (ply  drops  separated  by  20  ply  thicknesses)  with 

reinforced  resin  pockets. 


APPENDIX  B:  DIMENSIONAL  ANALYSIS 


B.l.  Introduction 

Dimensional  analysis  is  based  on  the  simple  idea  that  any  physical  law  is  independent 
of  the  arbitrary  chosen  units  of  measurement  used  for  describing  a specific  system.  Instead, 
any  physical  law  possesses  a fundamental  property  that  may  be  expressed  in  terms  of  variables 
intrinsic  to  the  system.  Dimensional  analysis  is  an  important  tool  in  developing  approximate 
response  functions  based  on  experimental  data  (either  physical  or  numerical  experiments). 
Dimensional  analysis  may  greatly  reduce  the  number  of  experiments  required  to  construct  the 
approximate  response  function,  by  reducing  the  number  of  variables  that  describe  the 
response.  The  variables  that  describe  the  response  are  known  as  governing  variables  (in  this 
dissertation  these  variables  are  also  referred  to  as  predictor  variables).  Dimensional  analysis  is 
discussed  in  detail  in  Barenblatt48  (1996,  pp.1-63),  which  makes  use  of  the  Buckingham  El- 
Theorem  to  formulate  a non-dimensional  response  quantity  in  terms  of  a group  of  non- 
dimensional  variables.  The  number  of  non-dimensionless  variables  is  in  general  smaller  than 
the  original  number  of  governing  variables.  The  Buckingham  El-Theorem  (Buckingham49 
1915)  may  be  stated  as: 

A physical  relationship  between  some  dimensional  (generally  speaking) 
quantity  and  several  dimensional  governing  variables  can  be  rewritten  as  a relationship 
between  some  dimensionless  quantity  and  several  dimensionless  products  of  the 
governing  variables;  the  number  of  dimensionless  products  is  equal  to  the  total 
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number  of  governing  variables  minus  the  number  of  governing  variables  with 
independent  dimensions 

An  independent  dimension  is  defined  as  a dimension  that  cannot  be  expressed  in  terms  of  a 
product  of  powers  of  the  remaining  variables. 

B.2.  Application  of  the  Buckingham  El-Theorem 
In  constructing  an  approximate  response  function  of  a system,  the  first  step  is  to 
express  the  response  in  terms  of  the  governing  variables  (predictor  variables)  as  follows: 

y = f{zl,z1,...,zk)  (B.l) 

The  governing  variables  may  then  be  divided  into  two  groups  such  that  the  m 
governing  variables  Z/  to  zm  have  independent  dimensions,  while  the  dimensions  of  the 
response  y and  remaining  (k  - m ) governing  variables  may  be  expressed  in  terms  of  the  m 
governing  variables  with  independent  dimensions,  as  follows: 

y = f(zl,...,zm,zm+l,...,zk)  (B.2) 

[y]  = [z,r,>[zj°<1)...[zj^) 

[zm+J  ] = kr™  [z2r™  ..\zm  m 

Note  that ye[  1,  (k  - m)\  and  that  the  independent  dimensions  account  for  all  the  fundamental 
dimensional  units  of  the  problem,  typically  force  [F],  length  [L]  and  possibly  time  [T]  for  a 
mechanics  problem. 

Based  on  Eq.  (B.3)  it  is  thus  possible  to  define  a non-dimensional  response  quantity  n 
and  a set  of  ( k-m ) non-dimensional  variables  n,  as  follows: 
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n = 


y 

z.amz*»...z  w 

i i m 


n.  =■ 


z “(m+JXl)  z °(»i+jX2)  2 “(m+jXm) 


(B.4) 


According  to  the  Buckingham  Id-Theorem  the  non-dimensional  response  quantity  is 
only  dependent  on  the  (k  - m)  non-dimensional  variables  as  follows: 

n = A(n  i , n 2 , . . . , n ) (B-5) 


B.3.  Example  Problem 

As  an  example  problem  consider  a cantilevered  beam  with  length  L and  circular  cross 
section  with  diameter  Das  shown  in  Fig  (B.l). 


L 


Figure  B.l:  Cantilevered  beam  example  problem. 

The  bending  stress  ctb  at  the  root  of  the  cantilevered  beam  may  then  be  written  in 
terms  of  the  governing  variables  (predictor  variables)  as: 

ctb  =f(P,L,D,y)  (B.6) 

where  P is  the  applied  load  and  y is  the  distance,  measured  from  the  neutral  axis  (see  Fig.  B.l), 
at  which  the  bending  stress  aB  is  evaluated.  The  fundamental  dimensional  units  of  the 
problem  are  force  [F]  and  length  [L].  The  applied  load  P (dimension  [F])  and  diameter  D 
(dimension  [L])  are  selected  as  the  governing  variables  with  independent  dimensions.  The 
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non-dimensional  bending  stress  aB  may  then  be  written  in  terms  of  two  non-dimensional 
variables  as  follows: 

^L=F(k  zl 

P \D'D) 

The  dimensional  analysis  thus  reduced  the  number  of  predictor  variables  from  four  to  only 
two.  In  terms  of  a quadratic  response  surface  approximation  the  reduction  in  predictor 
variables  would  reduce  the  number  of  parameters  in  the  full  response  surface  approximation 
from  15  to  only  6,  and  reduce  the  number  of  parameters  in  a quartic  response  surface 
approximation  from  35  to  only  10. 
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